Systems and methods for fast and precise frequency estimation

ABSTRACT

Systems and methods are provided for fast and precise estimation of frequency with relatively minimal sampling and relatively high tolerance to noise.

BACKGROUND

Frequency estimation systems are an important subsystem ofcommunication, navigation, radar and various other engineering systems.In some cases, the ability to perform relatively efficient and preciseestimation of frequency may be a critical component in system design.Indeed, an inability to perform relatively efficient and preciseestimation may significantly limit the performance of these systems asmeasured by various metrics of performance. For example, in high dynamicglobal positioning satellite (GPS) receiver applications, effectivelyacquiring and tracking the GPS carrier signal under dynamic conditionsmay limit the performance and applicability of these receivers tovarious applications. For this reason, various architectures suited tosuch applications have been previously proposed that acquire and trackthe GPS carrier under dynamic conditions, with the various differentmethods having a variety of limitations in terms of signal-to-noiseratio, initial frequency uncertainty, and other metrics of performance.

In terms of communication systems, precise frequency and phase may beimportant in communication systems involving coherent modulationtechniques such as multiple quadrature amplitude modulation (MQAM) andmultiple phase shift keying (MPSK). In certain communicationapplications, techniques, such as square law loop or Costas type loop,are used to derive the carrier frequency and phase from modulatedsignals, or a pilot signal is used which is tracked. The use of squarelaw loop or Costas type loop result in significant loss in terms ofphase noise of the reference carrier and phase ambiguity problemsassociated with phase ambiguity in the carrier phase equal to integermultiple of 2π/M for the MPSK signal. The use of pilot carrier resultsin a loss of signal power because a significant part of available poweris used up in the pilot. The ability to provide fast and accuratefrequency and phase estimates at very low signal-to-noise ratios (SNRs)may reduce the loss due to pilot carrier to an insignificant value.

More recently, precise and fast frequency acquisition and tracking havebecome increasingly important with the evolution of the OrthogonalFrequency Division Multiplexing (OFDM) in mobile communication systems.The OFDM modulation scheme may reduce problems of inter symbolinterference (ISI) caused by multipath propagation. It may also exhibitrelatively high performance in selective fading environments. Due tothese and other features, OFDM has become part of various standards suchas Worldwide Interoperability for Microwave Access (WiMax). Because OFDMis based on the orthogonality among various subcarrier signals, it isvery important that this orthogonality be maintained when thesesubcarriers are received at the receiver. However, the mobile wirelesschannels introduce frequency offsets which cause disruption of theorthogonality among the subcarriers resulting in mutual interferenceamong the various subcarriers. Therefore, it is desirable to preciselyestimate such frequency offsets and correct them to avoid problemsassociated with intercarrier interference. The offsets may be functionsof time and may vary with different subcarriers. Therefore, it isfurther desirable that precise estimation of the frequency offset bemade with a minimum requirement on the estimation time and SNR, which isalso limited in systems involving error correction coding techniques.

Some of the previous techniques for frequency estimation are based onthe extended Kalman filter (EKF). In these techniques, the state of theKalman filter is comprised of the signal phase, frequency and possiblyone or more derivatives of the frequency, and the measurement is anonlinear function of the state vector. In this approach, themeasurement function is linearized about the current estimate of thestate. Thus the technique based on EKF is more appropriate in thetracking mode when the initial estimate of the state is close to thetrue state and under a relatively high signal-to-noise ratio condition.However, when fast acquisition is required starting with a relativelyhigh frequency uncertainty and/or in low signal-to-noise ratioconditions, the EKF based methods may not meet the required performanceas they have relatively high thresholds for required signal power tonoise density ratio (P/N₀). Therefore, if the P/N₀ ratio is below thethreshold, the EKF may fail to converge and instead it may diverge. Inother words, when the EKF diverges, the state estimation error, insteadof converging to relatively small value as more and more measurementsamples are processed, increases with the number of samples processedand approaches a large estimation error. The EKF estimator also resultsin a phase locked loop (PLL) configuration with time-varying loop filtercoefficients and is thus also an “optimum” PLL. Thus the PLLs havesimilar limitations as the EKFs.

In fast Fourier transform (FFT) methods for frequency estimation, theN-point Fourier transform of a signal is evaluated, and the peak of theabsolute values of the FFT is searched. The FFT frequency correspondingto the peak is taken as the estimate of the unknown frequency associatedwith the signal. The FFT method has two limitations. First, theestimation error is of the order of the FFT resolution frequency that isequal to (2B/N) where 2B is the interval of frequency uncertainty.Second, in a relatively low SNR condition, selection of the peak among Nnoisy outputs of the transform involves detection errors correspondingto a large estimation error in the range of −B to B Hz. Increasing thesampling (N) reduces the resolution frequency; however, concurrently italso increases the probability of detecting the incorrect peak thusexhibiting threshold effects with respect to the SNR level.

The techniques that have been used in the estimation of the frequencyoffset in OFDM systems, based on the assumption that the offset is samefor all the sub-carriers, involve correlating the received signal with areference signal and thereby determining the relative phase between thetwo signals. By estimating the relative phase at two different timeinstances and dividing the difference between the two relative phasedifferences by the time difference provides the frequency estimate. Thecorrelation operation requires that the integration interval T_(I) forthe correlation operation may be much smaller than the inverse of theoffset frequency f_(a) else the signal amplitude at the correlatoroutput will be small and close to zero if T_(I)≅1/f_(a). On the otherhand, selection of a small T_(I) leads to relatively large noise at theintegrator output which may result in a noisy phase estimate.Differencing the noisy phase estimates may further accentuate theestimation error. Thus the correlation based approach has limitations interms of the magnitude of the frequency offset and the SNR conditionunder which it will provide a satisfactory result.

BRIEF DESCRIPTION OF THE DRAWINGS

The detailed description is set forth with reference to the accompanyingfigures. In the figures, the left-most digit(s) of a reference numberidentifies the figure in which the reference number first appears. Theuse of the same reference numbers in different figures indicates similaror identical items; however, various embodiments may utilize elementsand/or components other than those illustrated in the figures.

FIG. 1 is a block diagram illustrating an example architecture forproviding signal parameter estimation, in accordance with embodiments ofthe disclosure.

FIG. 2 is a block diagram illustrating an example architecture forproviding signal parameter estimation, in accordance with embodiments ofthe disclosure.

FIG. 3 shows a plot of example simulation results derived from thearchitectures of FIGS. 1 and 2, in accordance with embodiments of thedisclosure.

FIG. 4 shows a plot of example simulation results derived from thearchitectures of FIGS. 1 and 2, in accordance with embodiments of thedisclosure.

FIG. 5 shows a plot of example simulation results derived from thearchitectures of FIGS. 1 and 2, in accordance with embodiments of thedisclosure.

FIG. 6 shows a plot of example simulation results derived from thearchitectures of FIGS. 1 and 2, in accordance with embodiments of thedisclosure.

FIG. 7 shows a plot of example simulation results derived from thearchitectures of FIGS. 1 and 2, in accordance with embodiments of thedisclosure.

FIG. 8 shows a plot of example simulation results derived from thearchitectures of FIGS. 1 and 2, in accordance with embodiments of thedisclosure.

FIG. 9 shows a plot of example simulation results derived from thearchitectures of FIGS. 1 and 2, in accordance with embodiments of thedisclosure.

FIG. 10 shows a plot of example simulation results derived from thearchitectures of FIGS. 1 and 2, in accordance with embodiments of thedisclosure.

FIG. 11 shows a plot of example simulation results derived from thearchitectures of FIGS. 1 and 2 with comparisons to theoretical limits,in accordance with embodiments of the disclosure.

FIG. 12 shows a plot of example simulation results derived from thearchitectures of FIGS. 1 and 2 with comparisons to theoretical limits,in accordance with embodiments of the disclosure.

FIG. 13 shows a plot of example simulation results derived from thearchitectures of FIGS. 1 and 2 with comparisons to theoretical limits,in accordance with embodiments of the disclosure.

FIG. 14 shows a plot of example simulation results derived from thearchitectures of FIGS. 1 and 2 with comparisons to theoretical limits,in accordance with embodiments of the disclosure.

FIG. 15 shows a plot of example simulation results derived from thearchitectures of FIGS. 1 and 2 applied to a signal with a time-varyingfrequency, in accordance with embodiments of the disclosure.

FIG. 16 shows a plot of example simulation results derived from thearchitectures of FIGS. 1 and 2 applied to a signal with a time-varyingfrequency, in accordance with embodiments of the disclosure.

FIG. 17 shows a plot of example simulation results derived from thearchitectures of FIGS. 1 and 2 applied to a signal with a time-varyingfrequency, in accordance with embodiments of the disclosure.

FIG. 18 is a block diagram illustrating an example architecture forproviding signal parameter estimation, including frequency derivativeestimates and time-varying frequency estimates, in accordance withembodiments of the disclosure.

FIG. 19 shows a plot of example simulation results derived from thearchitecture of FIG. 18, in accordance with embodiments of thedisclosure.

FIG. 20 shows a plot of example simulation results derived from thearchitecture of FIG. 18, in accordance with embodiments of thedisclosure.

FIG. 21 is a block diagram illustrating an example architecture forproviding signal parameter estimation for multi-frequency signals, inaccordance with embodiments of the disclosure.

FIG. 22 is a block diagram illustrating an example architectureincluding a bank of filters for providing signal parameter estimationfor multi-frequency signals, in accordance with embodiments of thedisclosure.

FIG. 23 shows a plot of example simulation results derived from thearchitectures of FIGS. 21 and 22, in accordance with embodiments of thedisclosure.

FIG. 24 shows a plot of example simulation results derived from thearchitectures of FIGS. 21 and 22, in accordance with embodiments of thedisclosure.

FIG. 25 shows a plot of example simulation results derived from thearchitectures of FIGS. 21 and 22, in accordance with embodiments of thedisclosure.

FIG. 26 is a block diagram illustrating an example architecture forproviding signal parameter estimation using correlation detection, inaccordance with embodiments of the disclosure.

FIG. 27 is a block diagram illustrating an example second stage adaptiveKalman filter as included in the example architecture of FIG. 26, inaccordance with embodiments of the disclosure.

FIG. 28 is a block diagram illustrating an example correlation thresholddetector as included in the example architecture of FIG. 26, inaccordance with embodiments of the disclosure.

FIG. 29 shows a plot of example simulation results derived from thearchitecture of FIG. 26, in accordance with embodiments of thedisclosure.

FIG. 30 shows a plot of example simulation results derived from thearchitecture of FIG. 26, in accordance with embodiments of thedisclosure.

FIG. 31 shows a plot of example simulation results derived from thearchitecture of FIG. 26, in accordance with embodiments of thedisclosure.

FIG. 32 shows a plot of example simulation results derived from thearchitecture of FIG. 26, in accordance with embodiments of thedisclosure.

FIG. 33 shows a plot of example simulation results derived from thearchitecture of FIG. 26, in accordance with embodiments of thedisclosure.

FIG. 34 shows a plot of example simulation results derived from thearchitecture of FIG. 26, in accordance with embodiments of thedisclosure.

FIG. 35 shows a plot of example simulation results derived from thearchitecture of FIG. 26, in accordance with embodiments of thedisclosure.

FIG. 36 shows a plot of example simulation results derived from thearchitecture of FIG. 26, in accordance with embodiments of thedisclosure.

FIG. 37 shows a plot of example simulation results derived from thearchitecture of FIG. 26, in accordance with embodiments of thedisclosure.

FIG. 38 shows a plot of example simulation results derived from thearchitecture of FIG. 26, in accordance with embodiments of thedisclosure.

FIG. 39 is a block diagram illustrating an example architecture forproviding signal parameter estimation using correlation detection in alow signal-to-noise environment, in accordance with embodiments of thedisclosure.

FIG. 40 is a block diagram illustrating an example estimator, such as athird stage estimator of FIG. 39, in accordance with embodiments of thedisclosure.

FIG. 41 is a block diagram illustrating an example correlation block asincluded in the example architecture of FIG. 39, in accordance withembodiments of the disclosure.

FIG. 42 is a block diagram illustrating an example architecture forproviding signal parameter estimation using correlation detection in alow signal-to-noise environment for a multi-frequency signal, inaccordance with embodiments of the disclosure.

FIG. 43 shows a plot of example simulation results derived from thearchitecture of FIG. 39, in accordance with embodiments of thedisclosure.

FIG. 44 shows a plot of example simulation results derived from thearchitecture of FIG. 39, in accordance with embodiments of thedisclosure.

FIG. 45 shows a plot of example simulation results derived from thearchitecture of FIG. 39, in accordance with embodiments of thedisclosure.

FIG. 46 shows a plot of example simulation results derived from thearchitecture of FIG. 39, in accordance with embodiments of thedisclosure.

FIG. 47 shows a plot of example simulation results derived from thearchitecture of FIG. 39, in accordance with embodiments of thedisclosure.

FIG. 48 shows a plot of example simulation results derived from thearchitecture of FIG. 39, in accordance with embodiments of thedisclosure.

FIG. 49 shows a plot of example simulation results derived from thearchitecture of FIG. 39, in accordance with embodiments of thedisclosure.

FIG. 50 shows a plot of example simulation results derived from thearchitecture of FIG. 39, in accordance with embodiments of thedisclosure.

FIG. 51 shows a plot of example simulation results derived from thearchitecture of FIG. 39, in accordance with embodiments of thedisclosure.

FIG. 52 shows a plot of example simulation results derived from thearchitecture of FIG. 39, in accordance with embodiments of thedisclosure.

FIG. 53 illustrates an example computing system configured to executeany of the example architectures and/or algorithms disclosed herein, inaccordance with embodiments of the disclosure.

DETAILED DESCRIPTION

Embodiments of the disclosure are described more fully hereinafter withreference to the accompanying drawings, in which embodiments of thedisclosure are shown. This disclosure may, however, be embodied in manydifferent forms and should not be construed as limited to theembodiments set forth herein; rather, these embodiments are provided sothat this disclosure will be thorough and complete, and will fullyconvey the scope of the disclosure to those skilled in the art. Likenumbers refer to like elements throughout.

Embodiments of the disclosure are directed to, among other things,estimation of parameters associated with a signal. The disclosureprovides a variety of architectures based on Adaptive Kalman filter(AKF) and Fast Fourier Transform (FFT) methods. These variousarchitectures provide relatively fast and precise frequency estimatesalong with estimates of amplitudes, and phase for the case of single andmultiple carriers present in the received signal. The variousarchitectures operate on a vast range of signal-to-noise ratios whereinthose operating in low signal-to-noise ratio conditions involverelatively higher computational effort such that over a relatively largerange of SNRs, a high precision of estimation is achieved. The variousarchitectures presented also provide precise estimates for the caseinvolving significant frequency derivatives, wherein the frequencyderivatives may also be estimated from the received signal and accuratetracking of time-varying frequency is achieved.

The various architectures presented in this disclosure may further beused for single or multiple sinusoidal signals present in the sampledreceived signal. These architectures include one or more adaptive Kalmanfilters (AKF). The adaptive Kalman filter is comprised of two coupledKalman filters. The first Kalman filter (KF) is the state update Kalmanfilter (SUKF) that may estimate a state vector x(k) that may includein-phase and quadrature components of a sinusoidal signal received inthe presence of noise. The SUKF may provide the filtered estimates ofthe in-phase and quadrature components of the sinusoidal signal, thusfiltering out the additive noise from the measurements. The SUKF isbased on a state model, where the elements of a state transition matrixassociated with the state model are the cos and sin functions of thenormalized frequency ω₀=2πf_(a)/F_(s) where f_(a) is the frequency ofthe received signal and F_(s) is the sampling frequency. Based on theinformation provided by the SUKF, the second Kalman filter, termed theparameter update Kalman filter (PUKF), provides the estimates of aparameter vector β whose elements are cos (ω₀) and sin (ω₀). Based on adelayed filtered parameter vector estimate {circumflex over(β)}(k−1/k−1), the state transition matrix for the SUKF is generated.The two Kalman filters operate in a back to back manner providing thestate and the parameter vector estimates {circumflex over (x)}(k/k) and{circumflex over (β)}(k/k) for k=1, 2, . . . . From the state and theparameter vector estimates, the estimates for the frequency, phase, andamplitude are obtained in terms of simple functions of the state andparameter vector estimates. Unlike the EKF, the AKF does not involve anylinearization; therefore, it provides fast frequency estimates when theinitial frequency offset is relatively large and under medium SNRconditions. Simulations show that for an SNR of 7 dB, an RMS frequencyestimation error of about 3 Hz is obtained in about 10 samples when theinitial frequency is selected randomly from the range of −50 to +50 Hzand the sampling frequency is 400 Hz. With an increasing number ofsamples, the RMS error tends to be zero.

In the case of signals with randomly time-varying frequency withindependent incremental change at each sample, the AKF may closely trackthe frequency with the estimation error close to the random independentchange in frequency at each time instance. However, in the case ofsignificant frequency derivatives that may themselves be slowly varyingfunctions of time, the disclosure presents another embodiment whereinthe frequency derivatives are explicitly estimated, thereby achieving arapid acquisition followed by close tracking of frequency withrelatively small error. Another embodiment of the disclosure presentsacquisition and tracking of multiple, N sinusoidal signals present inthe received signal, wherein the various frequencies may also havesignificant frequency derivatives. In this embodiment, the estimatorreferred to as the adaptive Kalman filter-multiple frequency (AKF-MF),the state update Kalman filter (SUKF) may have dimension 2N and the PUKFmay have N KFs with each estimating a parameter vector of dimension 2.In yet another embodiment, the input signal is first filtered by a bankof M polyphase filters or the M FFT filters followed by the AKF-MFs thusresulting in the reduction of the dimension of the SUKF. The estimatesfrom the M AKF-MFs are appropriately combined to provide the finalestimates of the frequencies, frequency derivatives, amplitudes andphase of the N sinusoidal signals present in the received signal.

As the SNR is reduced below some value, there is a nonzero probabilitythat the estimates do not converge to small values that are of the orderof that predicted from the theoretical Cramer-Rao (C-R) bound. To avoidthis possibility, in another embodiment of the disclosure, the inputsignal, or a segment thereof, is correlated with the complex exponentialfunction of frequencies equal to the estimate provided by the AKF-MF. Ifthe correlation values are not above a certain bound in magnitude, theAKF-MF is rerun with different initial frequency estimates. The processis repeated until a threshold condition is satisfied or the number oftrials exceeds a specified number. At the end of the first estimationstep, the input signal is reprocessed with N individual single frequencyKFs, referred to herein as AKF-SFs, with their initial frequenciesprovided by the first stage estimation. The second stage may alsoinvolve processing of the input signal in segments and the operation offrequency shift at the end of processing each segment and reprocessingeach segment multiple times such that the frequency to be estimated bythe AKF-SF approaches zero at the end of each segment processing, andthe estimates from adaptive Kalman filter sub-systems (AKF-SSs) areaccumulated to provide final frequency estimates for the N inputfrequencies along with other parameters of interest. Simulation examplesshow that the use of this embodiment results in a root mean squared(RMS) frequency estimation error obtained with 50 random frequencyselections in the range of −100 to 100 Hz, and with 100 simulation runsfor each such frequency selection, of about 3 Hz for an SNR of 7 dB andin 10 input signal samples with no outliers. This compares to about a 2Hz RMS theoretical C-R lower bound on the RMS error.

Another embodiment of the disclosure, for the case of very low SNRs(˜−10 dB), involves three stages of estimation. In the first stage, theinput signal may be provided at the input to an M-point FFT with M≧Nwhere N is the number of frequencies in the input signal. The FFT mayprovide, at its output, an estimate of average power present at the MFFT frequencies. A frequency selector block may select N_(f)≦Mfrequencies from the M FFT frequencies corresponding to the highestaverage power levels. The selected N_(f) frequencies may be inputted toa bank of N_(f) AKF-SS that processes the input signal to arrive at thefinal frequency estimates f_(m) ^(K), m=0, 1, . . . , N_(f)−1. Theaverage power present at the N_(f) frequencies may be evaluated bycorrelation. For the case of N equal to 1, two frequencies out of theN_(f) frequencies f_(m) ^(K) with the highest power levels may beselected and their corresponding respective values may be compared tothreshold values. If the threshold condition is satisfied, the finalestimate f^(I) from the stage 2 of estimation is determined from the twoselected frequencies. If the threshold condition is not satisfied, thefirst and second processes of the estimation are repeated with afrequency offset introduced in the input signal. In the third process ofthe estimation, the input signal may be shifted in frequency by acomplex exponential function at frequency f^(I). The resulting frequencyshifted signal may be averaged over an appropriate interval, and thesignal thus obtained may be reprocessed by the AKF to achieve the finalestimate of the frequency f_(a). For the case of N>1, instead ofselecting two frequencies in the second stage, N frequencies with thehighest power may be selected. The N estimates from the third stage maybe provided to the third stage that is comprised of N AKFs, one for eachof the N frequencies. Simulations show that this architecture providesvery small frequency estimation error even at very low SNRs. Forexample, at an SNR of −10 dB, an RMS error equal to a fraction of 1 Hzis achieved in about 200 samples with the frequency selected randomlyfrom the range of −200 to 200 Hz. While the disclosure presents variouscomprehensive architectures for the estimation of frequency and otherparameters of multiple sinusoidal signals present in the input signal,these architectures can be adapted to a variety of applications thatinvolve frequency estimation in one or more forms.

Frequency Estimation with Adaptive Kalman Filter (AKF)

Various embodiments of the disclosure describe the adaptive Kalmanfilter (AKF) estimator and the FFT and correlation based estimators forfast and precise estimation of the frequencies, amplitudes and phases ofmultiple sinusoidal signals received in the presence of noise when theinterval of observation is relatively limited. The first embodiment ofthe disclosure describes the adaptive Kalman filter (AKF) estimator forthe case of a single sinusoidal signal present in the presence of noise.The frequency of the signal may remain constant over the interval ofobservation or may have some random time-varying component associatedwith it. In alternative embodiments of the AKF estimator, both thefrequency and the frequency derivative may be estimated when there is asignificant non-zero frequency derivative present during the observationinterval, and the derivative itself may be a constant or varying slowlyduring the interval of observation. Various embodiments of differentcomputational complexity are presented for signal-to-noise ratiosranging over medium to relatively very low values.

Estimation of Single Frequency

The embodiment of the disclosure for the AKF estimator for estimatingthe frequency, amplitude and phase of a signal, including a singlesinusoidal signal, is disclosed first. In this case, the intermediatefrequency signal of some unknown but constant frequency received in thepresence of band pass white Gaussian noise over an interval of time oflength T₀ is given byz _(IF)(t)=A cos (2π(f _(c) +f _(a))t+φ)+n(t);0≦t≦T ₀  (1)

In (1), f_(c) is known IF frequency, and the unknown frequency f_(a) hasa value in the range of −B to B Hz. The parameters A and φ are theamplitude and phase, respectively, both of which are also unknown. Theband pass white noise n(t) may have a one-side noise power spectraldensity of N₀ W/Hz. Thus the input signal to noise power ratio isA²/(4N₀B), and the signal power to noise power spectral density ratio isequal to A²/(2N₀) with the signal power A²/2. The corresponding baseband model obtained by the down conversion of the IF signal to complexbaseband is given byz _(i)(t)=A cos (ω_(a) t+φ)+n _(i)(t)  (2a)z _(q)(t)=A sin (ω_(a) t+φ)+n _(q)(t)  (2b)

In (2a, 2b) z_(i)(t) and z_(q)(t) are the in-phase and quadraturebaseband signals and ω_(a)=2πf_(a). The baseband noise processesn_(i)(t) and n_(q)(t) are independent white noise processes with thetwo-sided power spectral density N₀. This embodiment discloses anadaptive Kalman filter (AKF) estimator to estimate the unknownparameters A, ω_(a), and φ with minimum possible errors such that theperformance is close to theoretical bounds on such errors. Each of thebaseband signals has power A²/2 with the total baseband signal powerequal to A². In the embodiment, the signal power refers to the totalbaseband power.

The signals z_(i)(t) and z_(q)(t) are sampled at a sampling rate F_(s)that is at least two times the maximum possible value of the unknownfrequency f_(a)=ω_(a)/2π. The corresponding discrete time version of themeasurements is given byz _(i)(k)=A cos (kω ₀+φ)+n _(i)(k)  (3)z _(q)(k)=A sin (kω ₀+φ)+n _(q)(k);k=0,1, . . . ,N _(s)  (4)

In equations (3) and (4), k denotes the discrete time index, andω₀=ω_(a)T_(s) is the digital frequency in radians (rad), withT_(s)=1/F_(s), denoting the sampling period. The sampled noise sequencesn_(i)(k) and n_(q) (k) are independent random sequences with zero meanand variance σ_(n) ²=N₀/T_(s). For the purpose of estimation of theunknown parameters A, ω₀ and φ by the AKF, the state variables x₁(k) andx₂(k) are defined asx ₁(k)=A cos (kω ₀+φ)  (5)x ₂(k)=A sin (kω ₀+φ)  (6)

Using trigonometric expansion yields

$\begin{matrix}\begin{matrix}{{x_{1}\left( {k + 1} \right)} = {A\;{\cos\left( {{\omega_{0}\left( {k + 1} \right)} + \varphi} \right)}}} \\{= {{A\;{\cos\left( {{\omega_{0}k} + \varphi} \right)}{\cos\left( \omega_{0} \right)}} - {A\;\sin\;\left( {{\omega_{0}k} + \varphi} \right){\sin\left( \omega_{0} \right)}}}}\end{matrix} & (7)\end{matrix}$

From (5) (7), the following state update equation for the first statecomponent x₁(k) is obtainedx ₁(k+1)=β₁ x ₁(k)−β₂ x ₂(k)  (8)

where β₁ and β₂ are parameters given byβ₁=cos (ω₀);β₂=sin (ω₀)  (9)

Similarly expansion of x₂(k+1) yields

$\begin{matrix}\begin{matrix}{{x_{2}\left( {k + 1} \right)} = {A\;{\sin\left( {{\omega_{0}\left( {k + 1} \right)} + \varphi} \right)}}} \\{= {{A\;{\sin\left( {{\omega_{0}k} + \varphi} \right)}{\cos\left( \omega_{0} \right)}} + {A\;{\cos\left( {{\omega_{0}k} + \varphi} \right)}{\sin\left( \omega_{0} \right)}}}}\end{matrix} & (10)\end{matrix}$

From (5), (6), (8) and (10) one obtainsx ₂(k+1)=β₂ x ₁(k)+β₁ x ₂(k)  (11)

Equations (8) and (11) may be expressed in the following matrix vectorform

$\begin{matrix}{{x\left( {k + 1} \right)} = {{Fx}(k)}} & (12) \\{{{x(k)} = \begin{bmatrix}{x_{1}(k)} \\{x_{2}(k)}\end{bmatrix}};{F = \begin{bmatrix}\beta_{1} & {- \beta_{2}} \\\beta_{2} & \beta_{1}\end{bmatrix}}} & (13)\end{matrix}$

The measurement equations (3) and (4) for the AKF estimator may also bewritten in the matrix vector form

$\begin{matrix}{{{z(k)} = {\begin{bmatrix}{z_{i}(k)} \\{z_{q}(k)}\end{bmatrix} = {{{Hx}(k)} + {n(k)}}}};{{n(k)} = \begin{bmatrix}{n_{i}(k)} \\{n_{q}(k)}\end{bmatrix}}} & (14)\end{matrix}$

In (14), the measurement matrix H is equal to the 2×2 identity matrix.

FIG. 1 is a block diagram illustrating an example adaptive Kalman filter(AKF) 100 in accordance with embodiments of the disclosure. In oneaspect, the AKF 100 may be configured for the estimation of thefrequency f_(a), amplitude A and phase φ of the received signalz_(IF)(t) from the sampled noisy in-phase and quadrature signalsz_(i)(k) and z_(q)(k) in (3)-(4). The AKF may include two coupled Kalmanfilters. The first Kalman filter (KF) may be the state update Kalmanfilter (SUKF) 102 that may estimate the state vector x(k) that mayinclude the in-phase and quadrature components of the sinusoidal signalreceived in the presence of noise. The SUKF 102 may provide the filteredestimates of the in-phase and quadrature components of the sinusoidalsignal, thus filtering out the additive noise from the measurements. Theoutput of the SUKF 102 may be based on the state model, wherein theelements of the state transition matrix are the cos and sin functions ofthe normalized frequency ω₀=2πf_(a)/F_(s) where f_(a) is the frequencyof the received signal and F_(s) is the sampling frequency. Based on theinformation provided by the SUKF 102, the second Kalman filter, termedthe parameter update Kalman filter (PUKF) 106, may provide the estimatesof a parameter vector β with elements cos (ω₀) and sin (ω₀). In oneaspect, the SUKF 102 may provide its output, including, for example, thestate vector estimate {circumflex over (x)}(k/k) and the SUKF gain K(k)to the PUKF 106 via a state KF to parameter KF coupling 104. Based onthe delayed filtered parameter vector estimate {circumflex over(β)}(k−1/k−1), the state transition matrix for the SUKF 102 may begenerated.

The two Kalman filters 102, 106, therefore, operate in a back-to-backmanner providing the state and the parameter vector estimates{circumflex over (x)}(k/k) and {circumflex over (β)}(k/k) for k=1, 2, .. . . From the state and the parameter vector estimates, the estimatesfor the frequency, phase, and amplitude may be obtained in terms ofsimple functions of these estimates by an estimator of signal parameter110. Unlike an EKF, the AKF 100 does not involve any linearization.Therefore, it may provide fast frequency estimates when the initialfrequency offset is relatively large and under medium SNR conditions.

It will be appreciated that one or more of the constituent SUKF 102,PUKF 106, and other elements 104, 108, and 110 of the AKF 100 may beimplemented in a variety of mechanisms, including, for example, on aprocessor running instructions to implement the AKF 100. In one aspect,the elements of the AKF may be implemented as modules of a softwareapplication stored on one or more memories and accessible by one or moreprocessors configured to execute constituent instructions of the modulesand/or software application. Alternatively, each of the elements of theAKF may be implemented with discrete electrical elements communicativelycoupled therebetween. These elements may be digital elementsincorporating hardware, software, or both hardware and software.

The state update Kalman filter 102 equations for the estimation of thestate vector x(k) for the case of known matrix F are given by{circumflex over (x)}(k/k)={circumflex over (x)}(k/k−1)+K(k){tilde over(z)}(k/k−1)  (15a){tilde over (z)}(k/k−1)=z(k)−H(k){circumflex over (x)}(k/k−1)  (15b){circumflex over (x)}(k/k−1)=F(k−1){circumflex over (x)}(k−1/k−1)  (15c)S(k)=H(k)P(k/k−1)H ^(T)(k)+R(k)  (15d)K(k)=P(k/k−1)H ^(T)(k)S ⁻¹(k)  (15e)P(k/k)=P(k/k−1)−P(k/k−1)H ^(T)(k)S ⁻¹(k)H(k)P(k/k−1)  (15f)P(k+1/k)=F(k)P(k/k)F ^(T)(k)+G(k)Q(k)G ^(T)(k)  (15g)

In equations (15a-15g) {circumflex over (x)}(k/k) denotes the filteredestimate of the state x(k); {circumflex over (x)}(k/k−1) denotes thepredicted estimate of x(k) based on measurements up to time k−1, i.e.,z(0), z(1), . . . , z(k−1); {tilde over (z)}(k/k−1) is the predictionerror; and K(k) is the Kalman filter gain. The matrices P(k/k) andP(k/k−1) are the error covariance matrices for the filtered andpredicted estimates respectively of x(k). In general, the matrices F(k)and H(k) may be time-varying; however, in application to equations(12)-(14), these are independent of k, and the matrix G(k) in (15g) isequal to the identity matrix. The process covariance matrix Q is zerofor the equation (12). However, to increase the robustness of thealgorithm, Q may be set to some small non singular matrix. In anexponentially data weighted version of the Kalman filter, which has arobustness against certain model inaccuracies, the equation (15f) isreplaced byP(k+1/k)=λ⁻¹ F(k)P(k/k)F ^(T)(k)+G(k)Q(k)G ^(T)(k)  (15h)

In equation (15h), the parameter λ is the exponential data weightingcoefficient with 0<λ<1. The value of λ may be very close to 1, e.g.,0.99 and 1/(1−λ) is the approximate length of data memory. Thus withλ=0.99 the state update roughly depends only on the previous 100measurements.

From the definitions of the state vector x(k) in (5)-(6), the phase θ(k)at time k is given byθ(k)=tan₂ ⁻¹(x ₂(k),x ₁(k))  (16)

where tan₂ ⁻¹ denotes the 4 quadrant inverse tangent function, Similarlythe phase estimate is given in terms of the filtered estimate of thestate x(k) by{circumflex over (θ)}(k)=tan₂ ⁻¹({circumflex over (x)}₂(k/k),{circumflex over (x)} ₁(k/k))  (17)

The phase estimation error is obtained as{tilde over (θ)}(k)=θ(k)−{circumflex over (θ)}(k)  (18)

The filtered estimate of the amplitude A is also obtained from thefiltered estimate of the state x(k) asÂ(k)=[{circumflex over (x)} ₁ ²(k/k)+{circumflex over (x)} ₂²(k/k)]^(1/2)  (19)

A smoothed estimate of the amplitude is given by

$\begin{matrix}{{\hat{A}}_{S} = \left\lbrack {\frac{1}{N_{s} - N_{0}}{\sum\limits_{k = {N_{0} + 1}}^{N_{s}}{{\hat{A}}^{2}(k)}}} \right\rbrack^{1/2}} & (20)\end{matrix}$

where N_(s) is the total number of samples and N₀ is the number ofsamples that are discarded in the computation of the smoothed estimateof A.

The matrix F in the Kalman filter update equations (15) requiresknowledge of the parameters β₁ and β₂ which are unknown and need to beestimated. In the adaptive Kalman filter (AKF) 100 of the embodiment ofthis disclosure, these parameters are estimated from {circumflex over(x)}(k/k).

The equations (8) and (11) may be rewritten in the following equivalentform

$\begin{matrix}{{x_{1}\left( {k + 1} \right)} = {\left\lbrack {{x_{1}(k)} - {x_{2}(k)}} \right\rbrack\begin{bmatrix}\beta_{1} \\\beta_{2}\end{bmatrix}}} & (21) \\{{x_{2}\left( {k + 1} \right)} = {\left\lbrack {{x_{2}(k)}\mspace{14mu}{x_{1}(k)}} \right\rbrack\begin{bmatrix}\beta_{1} \\\beta_{2}\end{bmatrix}}} & (22)\end{matrix}$

Replacing the state x(k) by its filtered estimate {circumflex over(x)}(k/k), the equations (21) and (22) may be expressed in the followingmatrix vector form

$\begin{matrix}{{{{z_{p}(k)} \equiv \begin{bmatrix}{{\hat{x}}_{1}\left( {k/k} \right)} \\{{\hat{x}}_{2}\left( {k/k} \right)}\end{bmatrix}} = {{{H_{p}(k)}{\beta(k)}} + {n_{p}(k)}}};} & (23) \\{{H_{p}(k)} = \begin{bmatrix}{{\hat{x}}_{1}\left( {k - {1/k} - 1} \right)} & {- {{\hat{x}}_{2}\left( {k - {1/k} - 1} \right)}} \\{{\hat{x}}_{2}\left( {k - {1/k} - 1} \right)} & {{\hat{x}}_{1}\left( {k - {1/k} - 1} \right)}\end{bmatrix}} & (24)\end{matrix}$

In equations (23) and (24), the suffix on the variables z, H, n, etc.,signifies that these variables relate to the estimation of the parameterβ. The noise vector n_(p) accounts for the error in the estimate of thestate x(k) and has its covariance matrix equal to R_(p)(k) that isdetermined from the Kalman filter update equations (15). In case offixed frequency, the parameter vector β is not a function of k. However,to include a more general case, it may be assumed to vary according tothe following dynamics equationβ(k+1)=F _(p)β(k)+G _(p)(k)ξ_(p)(k)  (25)

where F_(p) and G_(p) are some appropriate matrices and ξ_(p)(k) denotesthe process noise, which may be zero mean and white, or otherwise spreadacross a relatively wide spectrum, with covariance matrix Q_(p). TheKalman filter update equations for the estimation of the parameter β(k)at the parameter update Kalman filter 106 from equations (23) and (25)is given by{circumflex over (β)}(k/k)={circumflex over (β)}(k/k−1)+K _(p)(k){tildeover (z)} _(p)(k/k−1)  (26a){tilde over (z)} _(p)(k/k−1)=z _(p)(k)−H _(p)(k){circumflex over(β)}(k/k−1)  (26b){circumflex over (β)}(k/k−1)=F _(p)(k−1){circumflex over(β)}(k−1/k−1)  (26c)S _(p)(k)=H _(p)(k)P _(p)(k/k−1)H _(p) ^(T)(k)+R _(p)(k)  (26d)K _(p)(k)=P _(p)(k/k−1)H _(p) ^(T)(k)S _(p) ⁻¹(k)  (26e)P _(p)(k/k)=P _(p)(k/k−1)−P _(p)(k/k−1)H _(p) ^(T)(k)S _(p) ⁻¹(k)H(k)P_(p)(k/k−1)  (26f)P _(p)(k+1/k)=F _(p)(k)P _(p)(k/k)F _(p) ^(T)(k)+G _(p)(k)Q _(p)(k)G_(p) ^(T)(k)  (26g)

In an exponentially data weighted version of the parameter update Kalmanfilter 106, the equation (26g) is replaced byP _(p)(k+1/k)=λ_(p) ⁻¹ F _(p)(k)P _(p)(k/k)F _(p) ^(T)(k)+G _(p)(k)Q_(p)(k)G _(p) ^(T)(k)  (26h)

In the Kalman filter equation (15) the matrix F(k−1) is replaced by itsestimate {circumflex over (F)}(k−1/k−1) given by

$\begin{matrix}{{\hat{F}\left( {k - {1/k} - 1} \right)} = \begin{bmatrix}{{\hat{\beta}}_{1}\left( {k - {1/k} - 1} \right)} & {- {{\hat{\beta}}_{2}\left( {k - {1/k} - 1} \right)}} \\{{\hat{\beta}}_{2}\left( {k - {1/k} - 1} \right)} & {{\hat{\beta}}_{1}\left( {k - {1/k} - 1} \right)}\end{bmatrix}} & \left( {26i} \right)\end{matrix}$

The two sets of equations (15) and (26) are updated in an interlacedmanner thus, the sequence of estimates are {circumflex over(β)}(k−1/k−1)→{circumflex over (F)}(k−1/k−1)→{circumflex over(x)}(k/k)→{circumflex over (β)}(k/k) for k=1, 2, . . . with some initialestimate for {circumflex over (β)}(0/0) based on some possibly zeroinitial estimate for the digital frequency ω₀.

In (26h), λ_(p) is the exponential data weighting coefficient with0<λ_(p)<1. The covariance matrix R_(p) in (26d) may be evaluated asfollows. The equations (21) and (22) may be combined to obtain thefollowing equation.

$\begin{matrix}{{{x(k)} = {{H_{pc}(k)}{\beta(k)}}};{{H_{pc}(k)} = \begin{bmatrix}{x_{1}\left( {k - 1} \right)} & {- {x_{2}\left( {k - 1} \right)}} \\{x_{2}\left( {k - 1} \right)} & {x_{1}\left( {k - 1} \right)}\end{bmatrix}}} & (27)\end{matrix}$

Subtraction of both sides of equation (23) from the respective sides ofequation (27) yieldsx(k)−{circumflex over (x)}(k/k)=[H _(pc)(k)−H _(p)(k)]β(k)−n_(p)(k);  (28)

From the definitions of the matrices H_(pc)(k) and H_(p)(k), theequation (28) may be written in the following alternative form

$\begin{matrix}{{{{\overset{\sim}{x}\left( {k/k} \right)} = {{\begin{bmatrix}{{\overset{\sim}{x}}_{1}\left( {k - {1/k} - 1} \right)} & {{\overset{\sim}{x}}_{2}\left( {k - {1/k} - 1} \right)} \\{{\overset{\sim}{x}}_{2}\left( {k - {1/k} - 1} \right)} & {{\overset{\sim}{x}}_{1}\left( {k - {1/k} - 1} \right)}\end{bmatrix}\begin{bmatrix}\beta_{1} \\\beta_{2}\end{bmatrix}} - {n_{p}(k)}}};}{\overset{\sim}{x}\left( {k/k} \right)} = {{x(k)} - {\hat{x}\left( {k/k} \right)}}} & (29)\end{matrix}$

With the definition of F in (13), one obtains from (29){tilde over (x)}(k/k)=F{tilde over (x)}(k−1/k−1)−n _(p)(k)  (30)Thusn _(p)(k)=F{tilde over (x)}(k−1/k−1)−{tilde over (x)}(k/k)  (31)

Now subtraction of (15a) from (12) yields{tilde over (x)}((k/k)=F{tilde over (x)}(k−1/k−1)−K(k){tilde over(z)}(k/k−1)  (32)

From (31) and (32) one obtains the covariance of n_(p)(k) asR _(p)(k)=E [n _(p)(k)n _(p) ^(T)(k)]=K(k)E [{tilde over(z)}(k/k−1){tilde over (z)} ^(T)(k/k−1)]K ^(T)(k)  (33)

Substitution for z(k) from (14) in (15b) results in{tilde over (z)}(k/k−1)=H(k){tilde over (x)}(k/k−1)+n(k)  (34)

As {tilde over (x)}(k/k−1) and n(k) are uncorrelated due to theindependence of noise sequence n(k), the covariance of the predictionerror {tilde over (z)}(k/k−1) is given byE[{tilde over (z)}(k/k−1){tilde over (z)} ^(T)(k/k−1)]=H(k)P(k/k−1)H^(T)(k)+R(k)=S(k)  (35)

Substitution of (35) in (33) yields the desired covariance matrix asR _(P)(k)=K(k)S(k)K ^(T)(k)=K(k)[H(k)P(k/k−1)H ^(T)(k)+R(k)]K^(T)(k)  (36)

For the case of F_(p)=I, corresponding to the case of time invariantparameter

and G_(p)=0, the Kalman filter equations (26) can be simplified toarrive at the weighted least squares algorithm. With the substitution ofF_(p)=I and G_(p)=0 and replacing k by (k−1) in (26h) yieldsP _(p)(k/k−1)=λ_(p) ⁻¹ P _(p)(k−1/k−1)  (37)

Substitution of (37) in (26f) yieldsP _(p)(k/k)=λ_(p) ⁻¹ P _(p)(k−1/k−1)−λ_(p) ⁻² P _(p)(k−1/k−1)H _(p)^(T)(k)×{λ_(p) ⁻¹ H _(p)(k)P _(p)(k−1/k−1)H _(p) ^(T)(k)+R _(p)}⁻¹ H_(p)(k)P _(p)(k−1/k−1)  (38)

The equation (38) can be organized in the following compact formP _(p)(k)=λ_(p) ⁻¹ [P _(p)(k−1)−P _(p)(k−1)H _(p) ^(T)(k)×{H _(p)(k)P_(p)(k−1)H _(p) ^(T)(k)+λ_(p) R _(p)}⁻¹ H _(p)(k)P _(p)(k−1)]  (39)

where in equation (39) Pp(k) is a more compact notation for P_(p)(k/k).The Kalman gain K_(p) for the case of constant β may be expressed asK _(p)(k)=P _(p)(k−1)H _(p) ^(T)(k)[H _(p)(k)P _(p)(k−1)H _(p)^(T)(k)+λ_(p) R _(p)]⁻¹  (40)

Multiplication of (39) by H_(p) ^(T)(k) results in

$\begin{matrix}{{{P_{p}(k)}{H_{p}^{T}(k)}} = {\lambda_{p}^{- 1}\begin{bmatrix}{{{P_{p}\left( {k - 1} \right)}{H_{p}^{T}(k)}} -} \\{{P_{p}\left( {k - 1} \right)}{H_{p}^{T}(k)}\left\{ {{{H_{p}(k)}{P_{p}\left( {k - 1} \right)}{H_{p}^{T}(k)}} + {\lambda_{p}R_{p}}} \right\}^{- 1} \times} \\{{H_{p}(k)}{P_{p}\left( {k - 1} \right)}{H_{p}^{T}(k)}}\end{bmatrix}}} & (41)\end{matrix}$

A simplification of (41) yieldsP _(p)(k)H _(p) ^(T)(k)=λ_(p) ⁻¹ P _(p)(k−1)H _(p) ^(T)(k){H _(p)(k)P_(p)(k−1)H _(p) ^(T)(k)+λ_(p) R _(p)}⁻¹ ×[H _(p)(k)P _(p)(k−1)H _(p)^(T)(k)+λ_(p) R _(p) −H _(p)(k)P _(p)(k−1)H _(p) ^(T)(k)]  (42)

From (40) and (42) one obtainsK _(p)(k)=P _(p)(k)H _(p) ^(T)(k)R _(p) ⁻¹  (43)

The weighted least squares algorithm is thus given byP _(p)(k)=λ_(p) ⁻¹ [P _(p)(k−1)−P _(p)(k−1)H _(p) ^(T)(k)×{H _(p)(k)P_(p)(k−1)H _(p) ^(T)(k)+λ_(p) R _(p)}⁻¹ H _(p)(k)P _(p)(k−1)]  (44a)K _(p)(k)=P _(p)(k)H _(p) ^(T)(k)R _(p) ⁻¹  (45b){circumflex over (β)}(k)={circumflex over (β)}(k−1)+K _(p)(k)[z_(p)(k)−H _(p)(k){circumflex over (β)}(k−1)]  (45c)

In (44), {circumflex over (β)}(k) is the compact notation for{circumflex over (β)}(k/k) in the case of constant β. From thedefinition of the parameters β₁ and β₂ in (9), the estimate of thedigital frequency ω₀ is given by{circumflex over (ω)}₀=tan₂ ⁻¹({circumflex over (β)}₂(k/k),{circumflexover (β)}₁(k/k))  (46)

with the estimation error in the estimate of the digital frequency ω₀given by{tilde over (ω)}₀(k)=ω₀−{circumflex over (ω)}₀(k)  (47)

The estimate and the estimation error for the analog frequency f_(a)given by{circumflex over (f)} _(a)={circumflex over (ω)}₀ F _(s)/(2π);{tildeover (f)} _(a) =f _(a) −{circumflex over (f)} _(a)  (48)

The set of equations (5), (6), (12)-(15), (17), (19) (20), (23)-(26),(33), (46) and (48) constitute the AKF algorithm for the estimation ofthe frequency, amplitude and phase of a single sinusoidal signalreceived in the presence of noise and is referred to as AKFS.

Referring now to FIG. 2, a relatively more detailed block diagramrepresentation of the example AKF 100 is shown. In this representation,the sampled signal z(k) defined in (14) is input to the state vectorupdate Kalman filter 102 that provides the filtered state estimate{circumflex over (x)}(k−1/k−1) at its output computed according toequations (15). The filtered state estimate {circumflex over(x)}(k−1/k−1) is pre-multiplied by the 2×2 matrix {circumflex over(F)}(k−1/k−1) to provide the predicted state estimate {circumflex over(x)}(k/k−1) according to (15c) wherein F(k−1) is replaced by itsestimate {circumflex over (F)}(k−1/k−1). The prediction error {tildeover (z)}(k/k−1) is evaluated according to (15b) by subtracting{circumflex over (x)}(k/k−1) from the measurement z(k) where H(k) in(15b) is substituted by the 2×2 identity matrix I₂. Multiplication of{tilde over (z)}(k/k−1) by the Kalman gain K(k) and adding the productto the predicted state estimate {circumflex over (x)}(k/k−1) results inthe filtered state estimate {circumflex over (x)}(k/k) at time k whichwhen delayed by one sample is equal to {circumflex over (x)}(k−1/k−1),thus completing the recursive loop for the update of the filtered stateestimate. The initial filter state {circumflex over (x)}(0/0) may beselected equal to [0 0]^(T).

Referring still to FIG. 2, the filter error covariance matrix P(k−1/k−1)is pre-multiplied by {circumflex over (F)}(k−1/k−1) and post-multipliedby the transpose of {circumflex over (F)}(k−1/k−1), and the resultingproduct is multiplied by 1/λ, where λ is the exponential data weightingcoefficient. The resultant product λ⁻¹{circumflex over(F)}(k−1/k−1)P(k−1/k−1){circumflex over (F)}^(T) (k−1/k−1) is added tothe matrix Q(k−1) according to (15h) where G(k)=I₂ and Q(k−1) isselected to be a constant diagonal matrix ε_(Q)I₂ for some smallpositive scalar ε_(Q). The innovation covariance matrix S(k) is obtainedby adding the result to the matrix R(k) that is equal to σ_(n) ²I₂ whereσ_(n) ² is the variance of the sampled noise sequences n_(i)(k) andn_(q) (k) in (14) and H(k)=I₂. Inverting the matrix S(k) and pre andpost multiplying the matrix S⁻¹(k) by P(k/k−1) and subtracting theresult from P(k/k−1) results in P(k/k) according to (15f). Thecovariance matrix P(0/0) is selected equal to εI₂ for some smallpositive scalar ε. Pre-multiplying S⁻¹(k) by P(k/k−1) yields the Kalmangain K(k) according to (15e).

The innovation covariance matrix S(k) is pre- and post-multiplied byK(k) and K^(T)(k) respectively providing the measurement covariancematrix R_(p)(k) for the parameter vector update Kalman filter accordingto (36). The elements of the vector {circumflex over (x)}(k−1/k−1)provide the measurement matrix H_(p)(k) given in (28), (29). The matrixH_(p)(k) is input to the parameter vector update Kalman filter for whichthe filtered state estimate {circumflex over (x)}(k/k) is themeasurement vector z_(p)(k).

The parameter vector update Kalman filter 106 is similar to the statevector update Kalman filter 102 with the filtered state estimate{circumflex over (x)}(k/k) replaced by the filtered parameter vectorestimate {circumflex over (β)}(k/k). For the case of constant parametervector β, the matrix F_(p)(k−1) is equal to I₂. The elements of thefiltered parameter vector estimate {circumflex over (β)}(k−1/k−1)provide the matrix {circumflex over (F)}(k−1/k−1) according to (26i)which is input to the state vector update Kalman filter. The two Kalmanfilters 102, 106 in FIG. 2 operate in a recursive and interlaced mannerproviding the estimates {circumflex over (β)}(k−1/k−1) and {circumflexover (x)}(k/k) for k=1, 2, . . . wherein the initial estimate{circumflex over (β)}(0/0) is set equal to [cos (ω⁰) sin (ω⁰)]^(T) withω⁰ denoting an a-priori estimate for the normalized frequency ω₀ whichmay be equal to 0.

From the filtered parameter vector estimate {circumflex over (β)}(k/k),the estimate for the normalized frequency ω₀ is obtained by the tan₂ ⁻¹() block that provides the four quadrant inverse tangent of the inputaccording to (46). The estimate for the analog frequency f_(a) may beobtained by multiplying {circumflex over (ω)}₀ by F_(s)/2π where F_(s)is the sampling frequency. The filtered state estimate {circumflex over(x)}(k/k) is input to a second tan₂ ⁻¹( ) block that yields the estimatefor the instantaneous phase θ(k) according to (17). The filtered stateestimate {circumflex over (x)}(k/k) is input to the norm ∥ ∥ block thatevaluates the Euclidean norm of the input providing the estimate Â(k)for the amplitude according to (19). From Â(k), the smoothed estimateÂ_(S) is obtained by the averaging block that obtains Â_(S) according to(20).

Simulation Results

As a non-limiting simulation example, consider the bandwidth B=200 Hzwith a sampling rate F_(s) selected equal to 400 Hz. The amplitude A=1volt with the bandpass noise one-sided noise spectral density N₀ isequal to 2.5×10⁻⁴ W/Hz with the signal power to noise spectral densityratio equal to 33 dB-Hz. This corresponds to the variance of thebandpass noise σ_(n) ²=2N₀B=0.1 W with the signal-to-noise ratio equalto 5 or 7 dB. The one-sided noise spectral density of both the basebandprocesses n_(i)(t) and n_(q)(t) is equal to 2N₀=5×10⁻⁴ W/Hz with thebaseband SNR for both the in-phase and quadrature components equal to 5(7 dB). The initial values for the matrices P, Q and P_(p) are selectedas

$\begin{matrix}{{{{P\left( {0/{- 1}} \right)} = \begin{bmatrix}5 & 0 \\0 & 5\end{bmatrix}};{{P_{p}\left( {0/{- 1}} \right)} = \begin{bmatrix}4 & 0 \\0 & 4\end{bmatrix}};}{{Q = \begin{bmatrix}0.001 & 0 \\0 & 0.001\end{bmatrix}};{{Qp} = 0};}} & (49)\end{matrix}$

The exponential data weighting coefficients are selected asλ=λ_(p)=0.999. FIGS. 3-10 illustrate example results obtained over theinterval of 25 samples. For the smoothed estimate of A, the number ofsamples N₀ is selected equal to 50.

FIG. 3 shows a plot of the two components of state estimation error{tilde over (x)}(k/k). FIG. 4 shows a plot of the phase estimation error{tilde over (θ)}(k) versus the number of samples k. FIG. 5 shows a plotof the error in the estimation of the analog frequency f_(a). The errorin the estimation of f_(a) may be relatively small after about 20samples, wherein the initial estimation error is 55 Hz. FIG. 6 shows aplot of the result ({tilde over (ω)}₀k)² as a function of k showing itsreduction with k. FIG. 7 shows a plot of the result for the trace (Rp)/2with Rp computed in (33). In the simulations, the matrix Rp has beenmodified by adding a diagonal matrix equal to 0.001 I₂ with I₂ denotingthe 2×2 identity matrix. The modified R_(p) matrix converges to 0.001I₂. FIGS. 8 and 9 show a plot of the results for an interval of 100samples showing that the errors in both the estimates of the state x(k)and the frequency f_(a) converge to zero. FIG. 10 shows a plot of thesmoothed estimate of the amplitude A from (19)-(20), with N₀ in (20)selected equal to 60, showing that the estimate is close to the truevalue of the amplitude equal to 1.0.

Comparison with Theoretical Bounds

The performance of the adaptive Kalman filter (AKF) 100 may be comparedto the theoretical bound on the frequency estimation error variance,which is also known as Cramer-Rao (C-R) bound. The Cramer Rao bound onthe variance of frequency estimation error {tilde over (f)}_(a) definedin (48) is given by

$\begin{matrix}{{\overset{\sim}{f}}_{{rm}\; s}^{CR} = {\left\{ {E\left\lbrack {\overset{\sim}{f}}_{a}^{2} \right\rbrack} \right\}^{1/2} \geq \left\{ {\frac{3}{2\pi^{2}}\left( \frac{2\sigma_{n}^{2}}{A^{2}} \right)\frac{F_{s}^{2}}{N_{s}^{3}}} \right\}^{1/2}}} & \left( {49a} \right)\end{matrix}$

The simulation examples compare the performance of the AKF algorithmwith the C-R bound. For the comparison, the AKF algorithm, as executedby the AKF 100 is simulated with the frequency f₀ selected randomly fromthe range (−B, B) Hz based on a uniform distribution and with theinitial frequency estimate equal to 0. For each such selection, the AKFalgorithm is simulated with n_(s)=100 independent noise sequences. Theprocedure is repeated with n_(f) different values of the frequency f₀selected on a random basis. For each selection of f₀, the RMS frequencyestimation error {tilde over (f)}_(rms) ^(i) as a function of the numberof samples processed is computed by the mean squared average over then_(s) simulation runs given by

$\begin{matrix}{{{{\overset{\sim}{f}}_{rms}^{i} = \left\{ {\frac{1}{n_{s}}{\sum\limits_{s = 1}^{n_{s}}\left\lbrack {\overset{\sim}{f}}_{i,s} \right\rbrack^{2}}} \right\}^{1/2}};{i = 1}},2,\ldots\mspace{14mu},n_{f}} & \left( {49\; b} \right)\end{matrix}$

FIG. 11 shows the RMS frequency estimation error as computed from (49b)for an initial estimation error of 55 Hz as a function of the number ofsamples processed. The initial rapid convergence of the RMS frequencyestimation error for one of the selections for f₀ shows that the RMSerror may be reduced to about 8 Hz from the initial error of 55 Hz inabout 8 samples. This is followed by a relatively slower convergencerate to an error of about 3 Hz in about 45 samples. FIG. 12 shows theRMS frequency estimation error as a function of the number of samplesfor the n_(f)=100 different frequency selections f₀. In FIG. 12, thesignal-to-noise ratio

$\gamma \equiv \frac{A^{2}}{2\;\sigma_{n}^{2}}$may be equal to 5 or 7 dB. The initial parameter error covariance matrixP_(p) is selected as

$\begin{matrix}{{P_{p}\left( {0/0} \right)} = \begin{bmatrix}1 & 0 \\0 & 1\end{bmatrix}} & \left( {49\; c} \right)\end{matrix}$

All the other parameters are selected, as in the non-limiting example,to compare with the C-R bound in (49a), and the RMS estimation errorevaluated in (49b) is averaged over all the n_(f)=100 frequencyselections as

$\begin{matrix}{{\overset{\sim}{f}}_{rms} = \left\{ {\frac{1}{n_{f}}{\sum\limits_{i = 1}^{n_{f}}\left\lbrack {\overset{\sim}{f}}_{rms}^{i} \right\rbrack^{2}}} \right\}^{1/2}} & \left( {49\; d} \right)\end{matrix}$

FIG. 13 shows a plot of the RMS frequency error obtained from the CramerRao bound in (49a) and that obtained from the AKF as evaluated in (49d).The AKF algorithm may provide an RMS error close to the theoreticalbound especially for a very low number of samples. For example, afterprocessing 4 samples, the AKF provides the RMS estimation error of 7.8Hz compared to the C-R bound of 6.2 Hz.

Time-Varying Frequency Estimation

In another non-limiting example, the frequency f_(a) is a randomlyvarying function of k according to the following equationf _(a)(k)=f _(a)(k−1)+δ_(f)(k)  (50)

The perturbation δ_(f)(k) is modeled as an independent zero meanGaussian sequence with variance σ_(f) ². From the definition of theparameters β₁ and β₂ in (9), the time-varying versions of theseparameters are approximated byβ₁(k)≅β₁(k−1)−2πT _(s)δ_(f)(k)β₂(k);β₂(k)≅β₂(k−1)+2πT_(s)δ_(f)(k)β₁(k)  (51)AsE[β ₁ ²(k)]=E[β ₂ ²(k)]=0.5  (52)

and δ_(f)(k) is independent of β₁(k−1) and β₂(k−1), one obtainsE[β ₁(k)−β₁(k−1)]² =E[β ₂(k)−β2(k−1)]²=0.5(2πT _(s))²σ_(f) ² ≡q_(a)  (53)

Thus the process covariance matrix Q_(p) in (26) is selected equal to

$\begin{matrix}{{Q_{p} = \begin{bmatrix}{q_{a} + ɛ} & 0 \\0 & {q_{a} + ɛ}\end{bmatrix}};{ɛ = {1.0 \times 10^{- 4}}}} & (54)\end{matrix}$

In the simulation example, σ_(f) ²=2 and Q is modified to 0.01 I₂ toaccount for added uncertainty in the estimate for F(k−1) required in theKalman filter equation (15). FIGS. 14-16 show the results for thisexample wherein all of the parameters are the same as for the case ofFIGS. 3-11. As may be observed from FIG. 14, the estimated frequencyremains close to the true frequency after an initial 10 samples. FIGS.15 and 16 plot the frequency estimation error and the two components ofthe state estimation error {tilde over (x)}(k/k) defined in (29)respectively. FIG. 17 plots the phase estimation error {tilde over(θ)}(k) defined in (18) as a function of k. After about 100 samples, thepeak-to-peak phase estimation error is about 0.4 rad, corresponding tothe RMS phase estimation error of about 0.07 rad. The RMS value of therandom change in frequency at any time k is equal to σ_(f)=√{square rootover (2)}Hz, and this component is independent of the measurements up totime k and therefore cannot be estimated at time k. Thus the idealestimation algorithm is expected to result in the RMS frequencyestimation error of σ_(f)=√{square root over (2)}Hz. From FIG. 15, thepeak frequency estimation error after the initial 10 samples is about 6Hz corresponding to an RMS error of about 2 Hz that is close to theideal case of √{square root over (2)} Hz with the difference of about0.6 Hz attributed to the estimation algorithm.

Estimation of Frequency Derivative

In another embodiment of the AKF estimator, both the frequency and thefrequency derivative may be estimated when there is a significantnon-zero frequency derivative present during the observation interval,and the derivative itself may be a constant or a slowly varying functionof time during the interval of observation.

In the case when the frequency is time-varying with its derivativeα_(a), the in-phase and quadrature baseband signals z_(i)(t) andz_(q)(t) are given by

$\begin{matrix}{{z_{i}(t)} = {{A\;{\cos\left( {\theta(t)} \right)}} + {n_{i}(t)}}} & (55) \\{{z_{q}(t)} = {{A\;{\sin\left( {\theta(t)} \right)}} + {n_{q}(t)}}} & (56) \\{{\theta(t)} = {{\omega_{a}t} + {\alpha_{a}\frac{t^{2}}{2}} + \varphi}} & (57)\end{matrix}$

From (57), the difference between the phase at time kT_(s) and(k−1)T_(s) is obtained as

$\begin{matrix}{{{\theta(k)} - {\theta\left( {k - 1} \right)}} = {{\omega_{a}T_{s}} + {\alpha_{a}\frac{t_{s}^{2}}{2}\left( {{2\; k} - 1} \right)}}} & (58)\end{matrix}$

With the time-varying normalized frequency ω₀(k) defined as

$\begin{matrix}{{{{\omega_{0}(k)} = {{\omega_{a}T_{s}} - {\alpha_{a}\frac{t_{s}^{2}}{2}} + {k\;\alpha_{a}T_{s}^{2}}}};{k = 1}},2,\ldots} & (59)\end{matrix}$

one obtainsθ(k)−θ(k−1)=ω₀(k)  (60)ω₀(k)=ω₀(k−1)+α₀;α₀≡α_(a) T _(s) ²  (61)

From (57), the instantaneous frequency ω_(i)(t) is given by

$\begin{matrix}{{\omega_{i}(t)} = {\frac{\mathbb{d}{\theta(t)}}{\mathbb{d}t} = {\omega_{a} + {\alpha_{a}t}}}} & (62)\end{matrix}$

From (59) and (62), the instantaneous frequency ω_(i)(t) at t=kT_(s) isrelated to the normalized frequency ω₀(k) byω_(i)(kT _(s))=[ω₀(k)+α₀/2]F _(s)  (63)

where F_(s) is the sampling frequency. In (61), α₀ may also be a slowlytime-varying function of k. The parameters β₁ and β₂ in (9) aretime-varying for this case and are given byβ₁(k)=cos (ω₀(k));β₂=sin (ω₀(k));k=1,2, . . .   (64)

With the state variables x₁(k) and x₂(k) defined asx ₁(k)=A cos (θ(k))  (65a)x ₂(k)=A sin (θ(k))  (65b)

the state space equations may be written as

$\begin{matrix}{{{x\left( {k + 1} \right)} = {{F(k)}{x(k)}}};{{x(k)} = \left\lbrack {{x_{1}(k)}{x_{2}(k)}} \right\rbrack^{T}}} & \left( {66\; a} \right) \\{{F(k)} = \begin{bmatrix}{\beta_{1}(k)} & {- {\beta_{2}(k)}} \\{\beta_{2}(k)} & {\beta_{1}(k)}\end{bmatrix}} & \left( {66\; b} \right)\end{matrix}$

The measurement equation is same as (14) rewritten here as (67).

$\begin{matrix}{{{z(k)} = {\begin{bmatrix}{z_{i}(k)} \\{z_{q}(k)}\end{bmatrix} = {{{Hx}(k)} + {n(k)}}}};{{n(k)} = \begin{bmatrix}{n_{i}(k)} \\{n_{q}(k)}\end{bmatrix}}} & (67)\end{matrix}$

In (67), H is equal to the 2×2 identity matrix. The state vector x(k) isestimated by the Kalman filter equations (15) with F(k−1) replaced byits estimate {circumflex over (F)}(k−1/k−1) obtained by replacing

₁k−

and β₂(k−1) in (66b) by their estimates {circumflex over (β)}₁(k−1/k−1)and {circumflex over (β)}₂ (k−1/k−1) respectively. The parameter vectorestimate {circumflex over (β)}(k/k−1)=[{circumflex over (β)}₁(k/k−1){circumflex over (β)}₂(k/k−1)]^(T) is obtained from the parameter updateequations (26) with z_(p)(k), H_(p)(k), and R_(p)(k) defined byequations (23), (24) and (36) respectively and with the matrix Q_(p) setequal to the 0 matrix, G_(p) equal to the (2×2) identity matrix and thematrix F_(p)(k−1) replaced by

$\begin{matrix}{{{\hat{F}}_{p}\left( {k - 1} \right)} = \begin{bmatrix}{\cos\left( {{\hat{\alpha}}_{0}\left( {k - 1} \right)} \right)} & {- {\sin\left( {{\hat{\alpha}}_{0}\left( {k - 1} \right)} \right)}} \\{\sin\left( {{\hat{\alpha}}_{0}\left( {k - 1} \right)} \right)} & {\cos\left( {{\hat{\alpha}}_{0}\left( {k - 1} \right)} \right)}\end{bmatrix}} & (68)\end{matrix}$

In (68), {circumflex over (α)}₀(k−1) denotes the estimate of α₀ based onthe measurement z(i) with i=1, 2, . . . , k−2. In order to obtain theestimate of {circumflex over (α)}₀ (k) for any k, the normalizedfrequency α₀(k) at time k is estimated by{circumflex over (ω)}₀(k)=tan₂ ⁻¹({circumflex over(β)}₂(k/k),{circumflex over (β)}₁(k/k));k=1,2 . . .   (69)

In (69), tan₂ ⁻¹( ) denotes the four quadrant inverse tangent functionsand {circumflex over (ω)}₀(0/0) may be set equal to an initial estimateω⁰ of ω(0) which may be equal to 0. In the parameter vector updateequation (26), {circumflex over (β)}(0/0) is correspondingly set equalto [cos (ω⁰)sin (ω⁰)]^(T). From (69), the first-order differences areobtained as{circumflex over (ω)}_(d)(k)={circumflex over (ω)}₀(k)−{circumflex over(ω)}₀(k−1);k=1,2, . . .   (70)

The estimates {circumflex over (ω)}_(d) (k) in (70) constitute the noisyestimate of the parameter α₀ wherein α₀ may also be a slowly varyingfunction of k. The smoothed estimate of α₀(k) is obtained by filteringthe estimates {circumflex over (ω)}_(d)(k) by an appropriate filter. Forexample, use of a first-order filter yields the following smoothedestimate of α₀(k).{circumflex over (α)}₀(k)=λ_(α){circumflex over(α)}₀(k−1)+(1−λ_(α)){circumflex over (ω)}_(d)(k);k=1,2, . . .   (71)

In (71), {circumflex over (α)}₀(0) may be set equal to an initialestimate of α₀ or may be set equal to 0. In (71), λ_(α) is a positiveconstant with 0<λ_(α)<1 and may determine the filter time constant whichis set in accordance with the expected rate of variation of α₀. Theparameter λ_(α) may itself be made a function of k; for example,initially λ_(α) may be set equal to 0.995 to provide a relatively smallweighting to the initial estimates {circumflex over (ω)}_(d)(k) with arelatively large estimation error and may be set to 0.99 after someinitial period to provide a fast tracking of α₀(k). The procedure can beextended to the case where a dynamic model for α₀(k) is specified interms of its derivatives by using a Kalman filter update equationsimilar to (26) and providing the estimate for α₀(k) and that of thederivatives of α₀(k) based on the specified dynamic model and treating{circumflex over (ω)}_(d)(k) as the noisy measurement of α₀(k).

Referring now to FIG. 18, an example architecture 1800 in accordancewith embodiments of the disclosure for the case of a relativelysignificant non-zero frequency derivative estimation along with thetime-varying frequency estimation is illustrated. The signal z(k) may beprovided to a frequency estimator block 1810 that may be similar to theAKF 100 of FIGS. 1 and 2. The frequency estimator block 1810 may includea state update Kalman filter 1812 and a parameter update Kalman filter1816. The frequency estimator block 1810 may provide at its outputs afiltered state estimate {circumflex over (x)}(k/k) and a filteredparameter vector estimate {circumflex over (β)}(k/k) from which afrequency estimate {circumflex over (f)}_(a)(k) may be obtained as inthe AKF 100 of FIGS. 1 and 2. Instead of being equal to the identitymatrix as may be the case for AKF 100 of FIGS. 1 and 2, may betime-varying and may be generated by a frequency derivative estimator1820.

The normalized frequency estimate {circumflex over (ω)}₀(k) may beprovided to a first-order differencing block 1830 which may generate adifference frequency {circumflex over (ω)}_(d)(k) according to (70). Thedifference frequency {circumflex over (ω)}_(d) (k) may be provided to afirst-order filter 1840 that generates at its output an estimate of thefrequency derivative {circumflex over (α)}₀(k) according to (71). Thedelayed version of {circumflex over (α)}₀ (k) may be input to thecascade of the cos ( ) and sin ( ) blocks and the scalar to vectorconverter 1850, followed by a matrix constructor block 1854 thatgenerates {circumflex over (F)}_(p)(k−1) from the vector [cos({circumflex over (α)}₀(k−1))sin ({circumflex over (α)}₀(k−1))]^(T)according to (68). The 2×2 matrix {circumflex over (F)}_(p)(k−1) may beinput to the parameter update Kalman filter 1816 for the generation ofparameter vector estimate {circumflex over (β)}(k/k).

A non-limiting simulation example is shown for the estimation of thetime-varying frequency with the frequency at time k=0 equal to 25 Hz andthe parameter α₀ equal to 0.007 corresponding to the parameter α_(a)equal to 0.007×400²=1120 rad/sec². This corresponds to a rate of changein frequency of about 178 Hz/sec. The initial estimates of the frequencyand the parameter α₀ at time k=0 are equal to −25 Hz and 0 respectively.The SNR for this example is 5 (7 dB) with the initial values for variousmatrices in the AKF, same as in (49). FIG. 19 shows a plot of thefrequency as a function of time along with its estimate obtained fromthe AKF. After about 100 samples, the frequency is tracked very closelyby its estimate. FIG. 20 plots the convergence of the parameter α₀showing a rapid initial convergence of the estimate to the trueparameter value. In (71) the filter parameter λ_(α) is initiallyselected equal to 0.995 and may be reduced to 0.99 after 150 samples.

In some embodiments of the disclosure, the first and higher orderderivatives may be obtained by a Kalman filter 3 that may replace thefrequency estimator block 1820 of FIG. 18. For example, for the case offirst and second derivative estimation, the state vector x_(h)(k) forthe Kalman filter 3 may be given byx _(h)(k)=[ω₀(k)α₀(k)β₀(k)]^(T)  (71a)wherein (71a) ω₀(k) and α₀(k) are the normalized frequency and thenormalized first derivative of frequency respectively and β₀=β_(a)T_(s)³ with β_(a) denotes the second derivative of the frequency ω(t). Thetime-variation of β₀ may be given byβ₀(k+1)=β₀(k)+n _(β)(k)  (71b)

In (71b), n_(β)(k) is a zero mean independent sequence with varianceσ_(β) ². The state transition equation for x_(h)(k) is given byx _(h)(k+1)=F _(h) x _(h)(k)+G _(h)ξ_(h)  (71c)where in (71c) G_(h) is the identity matrix and the covariance matrixQ_(h) of the noise vector ξ_(h) is given by

$\begin{matrix}{Q_{h} = {\begin{bmatrix}{1/20} & {1/8} & {1/6} \\{1/8} & {1/3} & {1/2} \\{1/6} & {1/2} & 1\end{bmatrix}\sigma_{\beta}^{2}}} & \left( {71\; d} \right)\end{matrix}$

The state transition matrix F_(h) in (71b) is given by

$\begin{matrix}{F_{h} = \begin{bmatrix}1 & 1 & {1/2} \\0 & 1 & 1 \\0 & 0 & 1\end{bmatrix}} & \left( {71\; e} \right)\end{matrix}$

The measurement equation for the Kalman filter 3 is given byz _(h)(k)=Hz _(h)(k)+n _(h)(k);H=[1 0 0]  (71f)with z_(h)(k) in (71f) made equal to the normalized frequency estimateobtained from the Kalman filter 2 state estimate {circumflex over(β)}(k/k) as in FIG. 18. Thus n_(h)(k)={tilde over (ω)}₀(k) and thecovariance of n_(h)(k) may be estimated from (51)-(53) and is given byR _(h)(k)=E[{tilde over (ω)} ₀ ²(k)]≅Trace[P _(p)(k/k)]  (71g)

On the basis of the state space model described by (71a-71g), the Kalmanfilter 3 may provide the estimate {circumflex over (x)}_(h) (k/k), thesecond component of which is the estimate of α₀(k) that may be inputtedto the block 1854 of FIG. 18 instead of that from block 1818 of FIG. 18.Derivatives of frequency of order higher than 2 can similarly beobtained by appropriate extension of the dimension of the state vectorx_(h)(k).

Estimation of Multiple Frequencies

In further embodiments of the disclosure, parameter estimation onmultiple sinusoidal signals received in the presence of noise may beperformed. In this case, the received intermediate frequency (bandpass)signal z_(IF)(t) is given by

$\begin{matrix}{{z_{IF}(t)} = {{\sum\limits_{m = 1}^{N}{A_{m}{\cos\left( {{2\;{\pi\left( {f_{c} + f_{am}} \right)}t} + \varphi_{m}} \right)}}} + {n(t)}}} & (72)\end{matrix}$

In (72), f_(IF) is the known intermediate frequency and the unknownfrequencies f_(am) have their values in the range of −B to B Hz for someB. The amplitudes A_(m) and phase φ_(m) associated with frequenciesf_(am) are also unknown. The bandpass white noise n(t) has a one-sidednoise power spectral density of N₀ W/Hz. The corresponding basebandmodel is given by

$\begin{matrix}{{z_{i}(t)} = {{\sum\limits_{m = 1}^{N}{A_{m}{\cos\left( {{\omega_{am}t} + \varphi_{m}} \right)}}} + {n_{i}(t)}}} & \left( {73\; a} \right) \\{{z_{q}(t)} = {{\sum\limits_{m = 1}^{N}{A_{m}{\sin\left( {{\omega_{am}t} + \varphi_{m}} \right)}}} + {n_{q}(t)}}} & \left( {73\; b} \right)\end{matrix}$

In (73a, 73b), z_(i)(t) and z_(q)(t) are the in-phase and quadraturesignals and ω_(arm)=2πf_(am) for m=1, 2, . . . , N. The basebandprocesses n_(i)(t) and n_(q)(t) are independent white noise processeswith the two-sided power spectral density N₀ W/Hz. It is required toestimate the parameters A_(m), ω_(am) and φ_(m) for m=1, 2, . . . , N.The signals z_(i)(t) and z_(q)(t) are sampled at a rate F_(s) that is atleast two times the maximum possible frequency B. The correspondingdiscrete time version of the measurements is given by

$\begin{matrix}{{z_{i}(k)} = {{\sum\limits_{m = 1}^{N}{A_{m}{\cos\left( {{k\;\omega_{m}} + \varphi_{m}} \right)}}} + {n_{i}(k)}}} & \left( {74\; a} \right) \\{{z_{q}(k)} = {{\sum\limits_{m = 1}^{N}{A_{m}{\sin\left( {{k\;\omega_{m}} + \varphi_{m}} \right)}}} + {n_{q}(k)}}} & \left( {74\; b} \right)\end{matrix}$

In equation (74), k denotes the discrete time index andω_(m)=ω_(am)T_(s), m=1, 2, . . . , N are the digital frequencies in radwith T_(s)=1/F_(s). The sampled noise sequences n_(i)(k) and n_(q)(k)are independent random sequences with zero mean and variance σ_(n)²=N₀/T_(s). In the case of multiple frequencies, the state vector x(k)is a 2N dimensional vectorx(k)=[x ₁(k)x ₂(k) . . . x _(2N)(k)]^(T)  (75)

With the components x_(m)(k) given byx _(2m-1)(k)=A _(m) cos (kω _(m)+φ_(m))x _(2m)(k)=A _(m) sin (kω _(m)+φ_(m));m=1, 2, . . . ,N  (76)

In the more general case where the m^(th) signal has a frequencyderivative α_(am), the state components in (76) are replaced byx _(2m-1)(k)=A _(m) cos (θ_(m)(k))  (77a)x _(2m)(k)=A _(m) sin (θ_(m)(k))  (77b)

where θ_(m)(k) is the sampled version of θ_(m)(t) given by

$\begin{matrix}{{\theta_{m}(t)} = {{\omega_{am}(t)} + {\alpha_{m}\frac{t^{2}}{2}} + \varphi_{m}}} & (78)\end{matrix}$

and as for the case of a single sinusoidal signal

$\begin{matrix}{{{\theta_{m}(k)} - {\theta_{m}\left( {k - 1} \right)}} = {\omega_{m}(k)}} & \left( {79\; a} \right) \\{{{{\omega_{m}(k)} = {{\omega_{am}T_{s}} - {\alpha_{am}\frac{T_{s}^{2}}{2}} + {k\;\alpha_{am}T_{s}^{2}}}};{k = 1}},2,\ldots} & \left( {79\; b} \right) \\{{{\omega_{m}(k)} = {{\omega_{m}\left( {k - 1} \right)} + \alpha_{m}}};{\alpha_{m} \equiv {\alpha_{am}T_{s}^{2}}}} & \left( {79\; c} \right)\end{matrix}$

The measurement equations (74a, 74b) may be written in the followingmatrix vector form

$\begin{matrix}{{{{z(k)} = {\begin{bmatrix}{z_{i}(k)} \\{z_{q}(k)}\end{bmatrix} = {{{Hx}(k)} + {n(k)}}}};}{{n(k)} = \begin{bmatrix}{n_{i}(k)} \\{n_{q}(k)}\end{bmatrix}}} & (80)\end{matrix}$

In (80), the measurement matrix H is a 2×2N matrix given in thefollowing partitioned formH[I ₂ I ₂ . . . I ₂]  (81)

where I₂ is the (2×2) identity matrix. The filtered and predicted stateestimates {circumflex over (x)}(k/k) and {circumflex over (x)}(k/k−1)are obtained from the Kalman filter equation (15) wherein thecorresponding covariance matrices P(k/k) and P(k/k−1) are squarematrices of dimension 2N. The process covariance matrix Q may be set tozero as for the case of single frequency estimation. The matrix F(k−1)in (15) is replaced by the estimate {circumflex over (F)}(k−1/k−1)obtained by the parameter vector update equations similar to (26). Thematrix {circumflex over (F)}(k−1/k−1) is a 2N×2N matrix given in thefollowing partitioned form

$\begin{matrix}{{\hat{F}\left( {k - {1/k} - 1} \right)} = {\quad\begin{bmatrix}{{\hat{F}}_{1}\left( {k - {1/k} - 1} \right)} & O & O & \ldots & O \\O & {{\hat{F}}_{2}\left( {k - {1/k} - 1} \right)} & O & \ldots & O \\\vdots & \vdots & \vdots & \ldots & \vdots \\O & O & O & \ldots & {{\hat{F}}_{N}\left( {k - {1/k} - 1} \right)}\end{bmatrix}}} & (82)\end{matrix}$

In (82), {circumflex over (F)}_(m)(k−1/k−1) for m=1, 2, . . . , N is a2×2 matrix given by

$\begin{matrix}{{{\hat{F}}_{m}\left( {k - {1/k} - 1} \right)} = \begin{bmatrix}{{\hat{\beta}}_{m\; 1}\left( {k - {1/k} - 1} \right)} & {- {{\hat{\beta}}_{m\; 2}\left( {k - {1/k} - 1} \right)}} \\{{\hat{\beta}}_{m\; 2}\left( {k - {1/k} - 1} \right)} & {{\hat{\beta}}_{m\; 1}\left( {k - {1/k} - 1} \right)}\end{bmatrix}} & (83)\end{matrix}$

The parameter vector estimate {circumflex over(β)}_(m)(k−1/k−1)=[{circumflex over (β)}_(m1)(k−1/k−1){circumflex over(β)}_(m2)(k−1/k−1)]^(T) is obtained recursively from the parameterupdate equation (26) with {circumflex over (β)}(k/k), {circumflex over(β)}(k/k−1), z_(p)(k), H_(p)(k) and R_(p)(k) in (26) replaced by{circumflex over (β)}_(m)(k/k), {circumflex over (β)}_(m)(k/k−1),z_(pm)(k), H_(pm)(k) and R_(pm)(k) respectively, where

$\begin{matrix}{{{z_{pm}(k)} = {\begin{bmatrix}{{\hat{x}}_{{2\; m} - 1}\left( {k/k} \right)} \\{{\hat{x}}_{2\; m}\left( {k/k} \right)}\end{bmatrix} = {{{H_{pm}(k)}{\beta_{m}(k)}} + {n_{pm}(k)}}}};} & (84) \\{{H_{pm}(k)} = \begin{bmatrix}{{\hat{x}}_{{2\; m} - 1}\left( {k - {1/k} - 1} \right)} & {- {{\hat{x}}_{2\; m}\left( {k - {1/k} - 1} \right)}} \\{{\hat{x}}_{2\; m}\left( {k - {1/k} - 1} \right)} & {{\hat{x}}_{{2\; m} - 1}\left( {k - {1/k} - 1} \right)}\end{bmatrix}} & (85)\end{matrix}$

The matrix Q_(p) is set equal to the zero matrix with G_(p) equal to the(2×2) identity matrix. The matrix R_(p)(k) in (33) is a 2N×2N matrix forthis case of multiple signals, with R_(pm) given by the following (2×2)submatrix of R_(p)(k).

$\begin{matrix}{{{R_{pm} = \begin{bmatrix}{R_{p}\left( {{{2\; m} - 1},{{2\; m} - 1}} \right)} & {R_{p}\left( {{{2\; m} - 1},{2\; m}} \right)} \\{R_{p}\left( {{2\; m},{{2\; m} - 1}} \right)} & {R_{p}\left( {{2\; m},{2\; m}} \right)}\end{bmatrix}};}{{m = 1},2,\ldots\mspace{14mu},N}} & (86)\end{matrix}$

The matrix F_(p)(k−1) in (26) is replaced by F_(pm)(k−1) given by

$\begin{matrix}{{{{F_{pm}\left( {k - 1} \right)} = \begin{bmatrix}{\cos\left( {{\hat{\alpha}}_{m}\left( {k - 1} \right)} \right)} & {- {\sin\left( {{\hat{\alpha}}_{m}\left( {k - 1} \right)} \right)}} \\{\sin\left( {{\hat{\alpha}}_{m}\left( {k - 1} \right)} \right)} & {\cos\left( {{\hat{\alpha}}_{m}\left( {k - 1} \right)} \right)}\end{bmatrix}};}{{m = 1},2,\ldots\mspace{14mu},N}} & (87)\end{matrix}$

In (87), {circumflex over (α)}_(m)(k−1) denotes the estimate of α_(m)based on the measurements z(i) for i=1, 2, . . . , k−2 and is estimatedas for the case of the single frequency estimation. Thus the normalizedfrequency ω_(m)(k) at time k is estimated by{circumflex over (ω)}_(m)(k)=tan₂ ⁻¹({circumflex over(β)}_(m2)(k/k),{circumflex over (β)}_(m1)(k/k));k=1,2 . . .   (88)

The normalized frequency estimate {circumflex over (ω)}_(m)(0/0) at time0 may be set to an initial estimate ω^(m) of ω_(m)(0) which may be equalto 0 and {circumflex over (β)}_(m)(0/0) in the parameter vector updateequation set equal to [cos (ω^(m)) sin (ω^(m))]^(T). From (88), thefirst order differences are obtained as{circumflex over (ω)}_(md)(k)={circumflex over (ω)}_(m)(k)−{circumflexover (ω)}_(m)(k−1);k=1,2, . . .   (89)

The estimate {circumflex over (ω)}_(md)(k) in (89) constitutes the noisyestimate of the parameter α_(m) which may also be a slowly varyingfunction of k. The smoothed estimate of ω_(m) may be obtained byfiltering the estimate {circumflex over (ω)}_(md)(k) by an appropriatefilter. For example, use of a first-order filter yields{circumflex over (α)}(k)=λ_(m){circumflex over(α)}_(m)(k−1)+(1−λ_(m)){circumflex over (ω)}_(md)(k);k=1, 2, . . .  (90)

In (90), {circumflex over (α)}_(m)(0) may be set equal to an initialestimate of α_(m) or may be set equal to 0. In (90), λ_(m) is a positiveconstant with 0<λ_(m)<1 and determines the filter time constant which isset in accordance with the expected rate of variation of α_(m). Theparameter, λ_(m) may itself be a function of k; for example, initiallyλ_(m) may be set equal to 0.999 to provide relatively small weighting tothe initial estimates {circumflex over (ω)}_(md)(k) with a relativelylarge estimation error and may be set to 0.99 after some initial periodto provide a fast tracking of α₀(k). The procedure may be extended tothe case where a dynamic model for α_(m)(k) is specified in terms of itsderivatives by using a Kalman filter update equation similar to (26) forproviding the estimate for α_(m)(k) and that of the derivatives of α₀(k)based on the specified dynamic model and treating {circumflex over(ω)}_(md)(k) as the noisy measurement of α₀(k).

From the definition of the components of the state vector in (77), thephase estimate at time k is obtained as{circumflex over (θ)}_(m)(k)=tan₂ ⁻¹({circumflex over (x)}_(2m)(k/k),{circumflex over (x)} _(2m-1)(k/k));m=1,2, . . . ,N  (91)

The filtered estimate of the amplitude A_(m) is given byÂ _(m)(k)=[{circumflex over (x)} _(2m-1) ²(k/k)+{circumflex over (x)}_(2m) ²(k/k)]^(1/2)  (92)

A smoothed estimate of the amplitude A_(m) may be obtained from (93).

$\begin{matrix}{{\hat{A}}_{mS} = \left\lbrack {\frac{1}{N_{s} - N_{0}}{\sum\limits_{k = {N_{0} + 1}}^{N_{s}}{{\hat{A}}_{m}^{2}(k)}}} \right\rbrack^{1/2}} & (93)\end{matrix}$

In (93), N_(s) is the total number of samples and N₀ is the number ofsamples that are discarded in the computation of the smoothed estimateof A_(m).

FIG. 21 illustrates another embodiment of the disclosure for theestimation of multiple frequencies present in an input signal z(k) withsignificant frequency derivatives using adaptive Kalman filtermulti-frequency (AKF-MF) 2100. The signal z(k) may be input to a stateupdate Kalman filter 2110 where the dimension of the state vector may beequal to 2N when N is the number of frequencies. The measurement matrixH at the input to the state update Kalman filter 2110 may be ofdimension 2×2N and is given by (81). The measurement noise covariancematrix R(k) at the input to the state update Kalman filter 2110 may begiven by σ_(n) ² I₂, and Q(k−1) may be equal to a constant diagonalmatrix ε_(Q)I_(2N) for a small positive scalar ε_(Q). The filtered stateestimate {circumflex over (x)}(k/k) may be provided to a vector splitter2120 that may partition the input vector {circumflex over (x)}(k/k) intoN vectors each of dimension 2 and denoted by {circumflex over(x)}¹(k/k), {circumflex over (x)}²(k/k), . . . {circumflex over(x)}^(N)(k/k) with{circumflex over (x)} ^(m)(k/k)=[{circumflex over (x)}_(2m-1)(k/k){circumflex over (x)} _(2m)(k/k)]^(T)  (93a)

A Kalman gain matrix K(k) of dimension 2N×2 may be input to the matrixsplitter to partition K(k) into N submatrices K¹(k), K²(k), . . . ,K^(N)(k) each of dimension 2×2 with

$\begin{matrix}{{{K^{m}(k)} = \begin{bmatrix}{K\left( {{{2\; m} - 1},1} \right)} & {K\left( {{{2\; m} - 1},2} \right)} \\{K\left( {{2\; m},1} \right)} & {K\left( {{2\; m},2} \right)}\end{bmatrix}}{{m = 1},2,\ldots\mspace{14mu},N}} & \left( {93\; b} \right)\end{matrix}$

In (93b), the argument k in the elements of the matrix has been droppedfor the notational simplification. The innovation covariance matrix S(k)may be pre- and post-multiplied by K^(m)(k) and its transpose K^(mT)(k),and the resulting product may be added to a constant matrix ε_(R)I₂, forsome small positive scalar ε_(R), providing a measurement covariancematrix R_(pm) for an m^(th) PUS (parameter update subsystem) block2130(A)-(N) (collectively referred to herein as 2130) for m=1, 2, . . ., N. The component state vector {circumflex over (x)}^(m) (k/k) may beprovided to the m^(th) PUS block 2130 to generate the estimates for thefrequency {circumflex over (f)}_(a,m) (k), frequency derivative{circumflex over (α)}_(a,m) (k), phase {circumflex over (θ)}_(m)(k) andamplitude Â_(m)(k) and the smoothed amplitude estimate Â_(m,S), for m=1,2, . . . , N. The operation of the PUS block 2130 has been explainedwith reference to block 1818 of FIG. 18 and, in the interest of brevity,will not be repeated here.

When the number of frequencies N in (72) is relatively large, thedimension of the covariance matrices P(k/k) and P(k/k−1) in the Kalmanfilter equation (15), equal to 2N×2N, is also correspondingly large. Inthis case, the matrix dimension can be reduced by filtering the complexvalued signal z_(c)(k)=z_(i)(k)+j z_(q)(k); j=√{square root over (−1)}by a set of M complex bandpass filters with their bandwidths equal to2B/M and their center frequencies selected so that the M filters havemutually disjoint pass bands. For example, with M=4, each of the 4filters has a bandwidth equal to B/2 with their respective centerfrequencies equal to −3B/4, −B/4, B/4, and 3B/4 respectively. Inalternative embodiments, there may be some non-zero overlap among theadjacent frequency bands of the M complex filters' output signals toavoid attenuating any frequency that may fall in a transition band ofthe filter.

The outputs of the M filters may be down converted to complex basebandso that each of the M down converted signals has its center frequencyequal to 0 and has frequencies in the range of −B/4 to B/4. The samplingrates of the down converted signals are decimated by a factor M to therate F_(s)/M. All of the operations of bandpass filtering, downconversion to complex baseband, and decimation can be performed by abank of polyphase filters.

The individual outputs of the polyphase filters are then processed bythe adaptive Kalman filters to acquire and track the frequencies presentin the respective outputs of the polyphase filters. The frequencyestimates thus obtained are adjusted by the center frequencies of thecorresponding bandpass filters to arrive at the frequency, frequencyderivative, and the amplitude estimates of the N signals present in thesignal z_(c)(k).

In an alternative embodiment, instead of using a polyphase filter, thefiltering of the signal z_(c)(i) may be performed by the FFT operationas follows. Let

$\begin{matrix}{{{{z^{m}(k)} = {\frac{1}{M}{\sum\limits_{i = {{({k - 1})}M}}^{{kM} - 1}{{z_{c}({\mathbb{i}})}{\mathbb{e}}^{{- {j2\pi}}\; f_{m}{\mathbb{i}}\; T_{s}}}}}};}{{f_{m} = {{- B} + {m\;\Delta\; f}}};}{{{\Delta\; f} = \frac{F_{s}}{M}};}{{m = 0},1,\ldots\mspace{14mu},{M - 1}}} & (94)\end{matrix}$

The operation in (94) constitutes frequency shifting of the signalz_(c)(i) by frequency f_(m) followed by a decimation by the factor Massumed to be even. The signal z^(m)(k) can also be obtained in terms ofthe discrete Fourier transform or the Fast Fourier transform (FFT) z_(f) ^(k) (of the k^(th) segment of z_(c)(i) of length M given by z^(k)=[z_(c)((k−1)M)z_(c)((k−1)M+1) . . . z_(c)(kM−1)] after beingmodified by multiplication of the i^(th) element of z ^(k) (by (−1)^(i),i.e., with

$\begin{matrix}{{{{{\overset{\_}{z}}_{f}^{k}(m)} = {\sum\limits_{i = 0}^{M - 1}{\left( {- 1} \right)^{i}{{\overset{\_}{z}}^{k}(i)}{\exp\left( {{- j}\; 2\;\pi\; i\;{m/M}} \right)}}}};}{{m = 0},1,\ldots\mspace{14mu},{M - 1}}} & (95) \\{{z_{c}^{m}(k)} = {\frac{1}{M}{{\overset{\_}{z}}_{f}^{k}(m)}}} & (96)\end{matrix}$

For a component z₀(k) at frequency f₀=(f_(m) ₀ +δf) contained in thesignal z_(c)(k) for some m₀, 0≦m₀≦M−1, with −Δf≦δf≦Δf, the k^(th)segment of that component complex valued signal is given byz ₀ ^(k)((i)=Ae ^(jφ)exp[j2π(f _(m) ₀ +δf)iT _(s) ];i=0,1, . . .M−1  (97)

where z ₀ ^(k)(i) denotes the i^(th) element of the segment z ₀ ^(k) andφ denotes the phase at the start of the k^(th) segment. From (94), thecorresponding output for the signal z ₀ ^(k) is given by

$\begin{matrix}{{z_{0}^{m}(k)} = {A\;{\mathbb{e}}^{j\;\varphi}{\sum\limits_{i = 0}^{M - 1}{\exp\left\lbrack {j\; 2\;{\pi\left( {f_{m_{0}} - f_{m} + {\delta\; f}} \right)}i\; T_{s}} \right\rbrack}}}} & (98)\end{matrix}$

From (98)

$\begin{matrix}{{{z_{0}^{m}(k)} = {\frac{A}{M}{\mathbb{e}}^{j\;\varphi}\exp\left\{ {j\;{\pi\left\lbrack {\frac{\delta\; f}{\Delta\; f}\frac{\left( {M - 1} \right)}{M}} \right\rbrack}} \right\}\frac{\sin\left( {{\pi\delta}\; f\text{/}\Delta\; f} \right)}{\left( {{\pi\delta}\; f\text{/}\left( {M\;\Delta\; f} \right)} \right)}}};{m = m_{0}}} & (99)\end{matrix}$

The expression for z₀ ^(m)(k) may also be approximated by

$\begin{matrix}{{z_{0}^{m}(k)} \cong \left\{ \begin{matrix}{{A\;{\mathbb{e}}^{j\;\varphi}\exp\left\{ {j\;{\pi\left\lbrack {\frac{\delta\; f}{\Delta\; f}\frac{\left( {M - 1} \right)}{M}} \right\rbrack}} \right\}\frac{\sin\left( {{\pi\delta}\; f\text{/}\Delta\; f} \right)}{\left( {{\pi\delta}\; f\text{/}\Delta\; f} \right)}};{m = m_{0}}} \\{0;{m \neq m_{0}}}\end{matrix} \right.} & (100)\end{matrix}$

As the adjacent frequencies f_(m) are spaced by Δf, for any frequency inthe interval f_(m) ₀ and f_(m) ₀ ₊₁ with f₀=f_(m) ₀ +δf; 0≦δf≦Δf, oneobtains

$\begin{matrix}{{{A_{1} \equiv {z_{0}^{m_{0}}(k)}} = {{Aa}_{1}{\mathbb{e}}^{j\;\varphi}{\mathbb{e}}^{j\;\psi_{1}}}};{a_{1} = \frac{\sin\left( {{\pi\delta}\; f\text{/}\Delta\; f} \right)}{\left( {{\pi\delta}\; f\text{/}\Delta\; f} \right)}};{\psi_{1} = {\left( {{\pi\delta}\; f\text{/}\Delta\; f} \right)\left( \frac{M - 1}{M} \right)}}} & \left( {101a} \right) \\{{{A_{2} \equiv {z_{0}^{m_{0} + 1}(k)}} = {{Aa}_{2}{\mathbb{e}}^{j\;\varphi}{\mathbb{e}}^{j\;\psi_{2}}}};{a_{2} = \frac{\sin\left( {\pi\left( {1 - {\delta\; f\text{/}\Delta\; f}} \right)} \right)}{\left( {\pi\left( {1 - {\delta\; f\text{/}\Delta\; f}} \right)} \right)}};{\psi_{2} = {\left( {\pi\left( {1 - {\delta\; f\text{/}\Delta\; f}} \right)} \right)\left( \frac{M - 1}{M} \right)}}} & \left( {101b} \right)\end{matrix}$

The adaptive Kalman filter for multiple frequencies (AKF-MF) 2100 may beapplied to signals z^(m)(k), m=0, 1, . . . , M−1. For any frequency f₀with f₀=f_(m) ₀ +δf; 0≦δf≦Δf, the AKF-MF 2100 applied to both z^(m) ⁰(k) and z^(m) ⁰ ⁺¹(k) result in frequency estimates f_(r1) and f_(r2)close to δf and −(Δf−δf) respectively with the estimate of δf given by

=(f_(r1)+f_(r2)+Δf)/2. From f_(r1) and f_(r2), the estimate of f₀ isgiven by f_(m) ₀ +(f_(r1)+f_(r2)+Δf)/2. The frequency estimates from anytwo AKF-MFs 2100 are combined in this manner if the estimates differ byless than the resolution frequency f_(res) which may be selected to be afraction of Δf. The corresponding estimates of the amplitudes Â₁ and Â₂provided by the two AKFM are expected to be close to the magnitudes ofA₁ and A₂ respectively given by (101).

Thus with the estimate of δf known, the coefficients a₁ and a₂ and thephase ψ₁ and ψ₂ may be evaluated from (101a, 101b) after replacing

by its estimate and the estimate of the amplitude A may be obtained bythe maximal ratio combining of Â₁ and Â₂ asÂ=|a ₁ Â ₁ e ^(−jψ) ¹ +a ₂ Â ₂ e ^(−jψ) ² |/(a ₁ ² +a ₂ ²)  (102a)

For the case of a2=0, the equation (102a) normalizes the amplitude A₁ bya₁ providing the correct estimate for the amplitude A. When the estimate

has significant error, an incoherent combining of the amplitudes Â₁ andÂ₂ may be used to yield the following amplitude estimate.Â=[a ₁ |Â ₁ +a ₂ |Â ₂|]/(a ₁ ² +a ₂ ²)  (102b)

The estimate for the phase is given by{circumflex over (φ)}=[arg(A ₁)−ψ₁+arg(A ₂)−ψ₂]/2  (103)

FIG. 22 illustrates a further example embodiment of the disclosurewherein the input signal, including multiple frequencies, may be splitinto a number M of signals that have approximately disjoint spectrumwith some possible overlap using the architecture 2200 of FIG. 22. Thecomplex valued signal z_(c)(i), where i denotes input time index, may beprovided to a bank of FFT filters 2210. In particular, z_(c)(i) may beinput to a serial to parallel converter 2214 that may split the inputsignal z_(c)(i) into M parallel streams denoted by z ^(k)(0), . . . z^(k)(M−1). The signals with odd index in their arguments may beinverted. For any time index k, the signals after being multiplied by ±1may be input to an FFT block 2222 that may provide the FFT of the inputvector z _(f) ^(k)(m), m=0, 1, . . . , M−1 as the output of the FFTblock 2222 according to (95). The FFT outputs after normalization by Mmay include filtered signals z_(c) ^(m)(k), m=0, 1, . . . , M−1; k=0, 1,. . . defined in (94). The signal z_(c) ^(m)(k) may be input to acomplex to real converter 2230(A)-(N) (collectively referred to hereinas 2230) that outputs a dimension 2 vector z^(m)(k) with its elementsequal to the real and imaginary parts, respectively, of the signal z_(c)^(m)(k).

The two dimensional vector z^(m)(k) is input to the AKF-MF block that isalso inputted with the matrices H_(m)(k), Q_(m)(k−1), and R_(m)(k). Thematrix H_(m)(k) is given by (81) and is of dimension 2×2N_(m),Q_(m)(k−1)=ε_(Q)I_(2Nm) for some small positive scalar ε_(Q) andR_(m)(k)=σ²I₂ with σ²=σ_(n) ²/M and with N_(m) denoting the number offrequencies in the signal z_(c) ^(m)(k) and may be replaced by someappropriate upper bound. Partitioning the input signal z_(c)(i) in thismanner reduces both the dimension of the state vector and various othermatrices and at the same time reduces the noise variance by the factorM.

Referring still to FIG. 22, the AKF-MF block 2250(A)-(N) (collectivelyreferred to herein as 2250) may provide, as its output, the frequencyestimates {circumflex over (f)}_(m,i)(k), frequency derivative{circumflex over (α)}_(m,i)(k) (not shown), phase {circumflex over(φ)}_(m,i)(k) and the amplitude Â_(m,i)(k), for i=1, 2, . . . , N_(m)for m=0, 1, . . . , M−1. The operation of the AKF-MF block 2250 has beenexplained with reference to FIG. 21, and in the interest of brevity,will not be repeated here. The estimates of frequencies after beingmodified by the frequencies (−B+mΔf), frequency derivatives, amplitudes,and phase may be input to a collator/combiner block 2270. In thecollator/combiner block 2270, any two frequency estimates in any twoadjacent AKF-MF blocks 2250 that are within the resolution frequency ofeach other are combined into one estimate obtained by either a simplearithmetic average of the two frequencies or by a weighted average withthe weights determined by the corresponding amplitudes with similaraveraging of the frequency derivatives. The corresponding amplitudes andphase may be combined according to (102)-(103). For the case of constantfrequencies, the proximity may be determined on the basis of only thefinal estimates at time N_(s) where as for the time-varying frequencies,the two frequency estimates should be close to each other over aninterval of time to be combined into one. The collator/combiner block,after possibly combining frequency estimates, may provide the Nfrequencies, frequency derivatives, amplitudes and phase as the output.In an alternative embodiment, the operation of splitting the inputcomplex signal z_(c)(k) into M parallel streams denoted by z _(c) ⁰(k),z ¹(k), . . . z ^(M-1)(k) with approximately mutually disjointed spectramay be performed by a bank of polyphase filters instead of the bank ofFFT filters 2210 as shown in FIG. 22.

Simulation Results

FIG. 23 shows a plot of the frequency estimation error for the case ofN=2 frequencies at 76 Hz and 135 Hz with their respective amplitudes A₁and A₂ both equal to 0.5. The noise variance σ_(n) ²=0.1 corresponds tothe total signal power to noise power ratio of 2.5 (3.5 dB). The matrixP(0/0) is equal to 5I₄ with I₄ denoting the 4×4 identity matrix andP_(p)(0/0)=0.1I₂, Q=1.0×10⁻⁴I₄. The initial estimates for the twofrequencies are 55 Hz and 161 Hz, respectively. The estimation error inthe estimates of both the frequencies may become very small with about40 samples. FIG. 24 shows a plot of the error in the amplitude estimatesÂ₁ and Â₂ obtained from (92). FIG. 25 shows a plot of the amplitudeestimation errors for the case when both the amplitudes A₁ and A₂ areequal to 1 corresponding to the total signal power to noise power ratioof 10 (10 dB). In this case, the amplitude estimation errors are alsosmall. The frequency estimation errors for this case are relativelysmall. The errors in the smoothed estimates of amplitudes Â_(1S) andÂ_(2S) obtained from (93) may be even smaller than those in FIGS. 24 and25.

Adaptive Kalman Filter-Correlation Detector (AKF-CD)

When the signal-to-noise ratio is relatively small and the initialfrequency uncertainty is relatively large, there may be some smallprobability that the AKF estimation algorithm may not converge to thetrue frequency with relatively small estimation error. To avoid such apossibility, a correlation detector is used to validate the frequencyestimates provided by the AKF.

In the AKF-CD architecture, the estimation interval of N_(s) samples issegmented into m_(s) segments of n_(s) samples each wherein in thespecial case n_(s) may be equal to N_(s). In this architecture, the AKFalgorithm with various matrices described in equations (74) to (103) maybe run for n_(s) samples to estimate {circumflex over (ω)}_(m) (k) and{circumflex over (α)}_(m) (k) for k=1, 2, . . . , n_(s); m=1, 2, . . . ,N where N is the number of frequencies. From (79), the phase estimate{circumflex over (θ)}_(m)(k) is obtained as{circumflex over (θ)}_(m)(k)={circumflex over (θ)}_(m)(k−1)+{circumflexover (ω)}_(m)(k);k=1,2, . . . ,n _(s)  (104)

In (104), {circumflex over (θ)}_(m)(0) may be set to 0. The complexvalued signal z_(c)(k)=z_(i)(k)+jz_(q)(k); j=√{square root over (−1)} iscorrelated with the signal s_(m)(k) in a semi-coherent manner, wheres _(m)(k)=exp(−j2π{circumflex over (θ)}_(m)(k));k=1,2, . . . ,n_(s)  (105)

For the semi-coherent correlation, the interval n_(s) is divided intom_(c) segments of length n_(c), and the coherent correlation is computedfor each of the m_(c) intervals as

$\begin{matrix}{{{{r_{c}^{m}(l)} = {\sum\limits_{k = {{{({l - 1})}n_{c}} + 1}}^{l\mspace{14mu} n_{c}}\;{{z_{c}(k)}{s_{m}(k)}}}};{l = 1}},2,\ldots\mspace{14mu},m_{c}} & (106)\end{matrix}$

In (106), the interval n_(c) is determined such that some bound on theexpected frequency estimation error in Hz obtained after processing then_(s) samples is much smaller than the F_(s)/n_(c) Hz wherein the boundon the expected frequency estimation error in Hz, for example, may besome multiple of the C-R bound. From (106), the correlation over theinterval of n_(s) samples is obtained as

$\begin{matrix}{r_{m}^{i} = \left\{ {\sum\limits_{l = {l_{0} + 1}}^{m_{c}}\;\left| {r_{m}^{c}(l)} \right|^{2}} \right\}^{1\text{/}2}} & (107)\end{matrix}$

In (107), l₀ that may be 0 denotes the number of initial n_(c) lengthsegments that are discarded in the computation of the incoherentcorrelation r_(m) ^(i). The number of samples l₀n_(c) accounts for theinitial transient period where the estimation error is relatively high.For the case of α_(m)=0, where the frequency ω(k) is constant, thecoherent correlation in (106) reduces to

$\begin{matrix}{{{{r_{m}^{c}(l)} = {\sum\limits_{k = {{{({l - 1})}n_{c}} + 1}}^{l\mspace{14mu} n_{c}}\;{{z_{c}(k)}{\exp\left( {{- j}\; k\mspace{14mu}{\hat{\omega}}_{m}} \right)}}}};{l = 1}},2,\ldots\mspace{14mu},m_{c}} & (108)\end{matrix}$

In (108), {circumflex over (ω)}_(m) denotes the estimate of thenormalized frequency in radians based on the processing of n_(s)samples. In this case, l₀ in (107) may be set to 0. In the ideal case ofsmall frequency error and relatively small noise variance, thecorrelation in (107) is nearly equal to r_(m) given byr _(m) ={A _(m) ²+2σ_(n) ² /n _(c)}^(1/2) ≅A _(m)  (109)

In the correlation detector, r_(m) ^(i) computed from (107) is comparedto a correlation threshold γ_(m) ^(th) which may be equal to somefraction of an estimate of A_(m). If r_(m) ^(i) is greater than or equalto the threshold, the estimate is deemed to converge to some smallerror. If all the N frequency estimates do not satisfy the thresholdconditions, the AKF algorithm described in equations (74) to (103) isrepeated for the same initial n_(s) samples but with different randomfrequency initializations, wherein the initializations for thosefrequencies which satisfy the threshold conditions are set equal totheir estimates obtained in the first run. The frequencies for which thethreshold condition is not satisfied are initialized randomly over someappropriate interval; for example, over the interval (−B/4, B/4). In theabsence of estimates of the individual A_(m), the threshold may beapplied to the total power estimated as

$\begin{matrix}{P_{T} = {\sum\limits_{m = 1}^{N}\;\left( r_{n}^{i} \right)^{2}}} & (110)\end{matrix}$

The power estimated in (110) is compared with some fraction of its priorestimate P_(th) which may, for example, be obtained as

$\begin{matrix}\left. {P_{th} \cong {\frac{1}{n_{s}}\sum\limits_{i = 1}^{n_{s}}}}\; \middle| z_{c} \middle| {}_{2}{{- 2}\;\sigma_{n}^{2}} \right. & (111)\end{matrix}$

If P_(T)>δP_(th) for some selected value of δ with 0<δ≦1, then the firstpart of the estimation procedure may be complete. Otherwise, the AKFalgorithm in equations (74) to (103) may be repeated for the sameinitial n_(s) samples but with different random frequencyinitializations until either the threshold condition is satisfied or thenumber of trials exceed some specified threshold N_(t). A higher valueof δ may in general require a higher number of trials.

Once the threshold condition is satisfied, each of the N frequencies istracked by further processing of the measurements. In this phase of thealgorithm, the second set of n_(s) complex measurements is shifted infrequency by the estimate {circumflex over (ω)}_(m)(n_(s)) for them^(th) normalized frequency ω_(m)(n_(s)) obtained at the end of thefirst phase according to equation (112).z _(c,m)(k)=z _(c)(k)exp(−jk{circumflex over (ω)} _(m)(n _(s)));k=n_(s)+1, . . . ,2n _(s)  (112)

The AKF algorithm processes the n_(s) samples of the measurementsz_(cm)(k) to obtain the estimates of the frequency (ω_(m)(k)−{circumflexover (ω)}_(m)(n_(s)), k=n_(s)+1, . . . , 2n_(s). In processing thesecond set of n_(s) measurements, the parameter error covariance matrixP_(p)(k/k−1) may be increased by an appropriate value so as to increasethe convergence rate. The estimates of the frequencies (ω_(m)(k)−{circumflex over (ω)}_(m)(n_(s)) thus obtained are adjusted by{circumflex over (ω)}_(m)(n_(s)) to account for the frequency shiftproviding the estimates for the frequency {circumflex over (ω)}_(m)(k),k=n_(s)+1, . . . , 2n_(s). This procedure of shifting the frequency ofthe complex signal by the frequency estimate obtained at the end of theprevious segment and applying the AKF to process the current segment andmodifying the resulting estimates by the amount of frequency shift isrepeated for the remaining m_(s)-2 segments. A further improvement inthe estimates may be obtained by reprocessing all of the m_(s) segmentsN_(r)-1 times in a manner similar to the processing of the m_(s)-1segments in the first run where N_(r) denotes the total number of runswith N_(r)≧1. For the case of time-varying frequency, the amount offrequency shift at the start of the second and consecutive runs needs tobe appropriately adjusted.

In case of multiple frequencies, the signal z_(c)(k) may be segmentedafter the initial phase into M signals by the use of a polyphase filteror using FFT as in (94)-(96) such that the number of frequencies in anyone segment is small and the dimension of the state vector x(k) in theAKF-MF is relatively small for processing any one segment.

FIG. 26 is a block diagram illustrating another embodiment of thedisclosure architecture 2600, where a correlation detection is performedto assess if the estimated frequency is close to the true frequency asdescribed in (104)-(111). If the correlation threshold condition is notsatisfied, the input signal may be processed again starting with a newinitial estimate for the frequency. The input signal may be stored in amemory 2610 so that if required the measurements can be reprocessed.When the number of available samples is not highly restricted, thestorage of the input samples is not necessary, because the iterationscan be performed with real time input samples. The architecture 2600 mayinclude a switch S₁, depending upon the state of the control signal I₆,that permits processing of either the real time input data or that fromthe memory. Switch S₂, depending upon the state of the control signalI₇, may route the complex signal z_(c)(k) to either a first stage AKF-MF2630 via a complex to real converter 2620 or to second stage processing.Switch S₃, depending upon the state of the control signal I₈, may routethe complex signal z_(c)(k) to either a bank of correlation detectorsCTD1, . . . , CTDN 2640(1)-(N) (collectively referred to herein as 2640)or to the bank of stage 2 adaptive Kalman filters AKF1, . . . , AKFN2650(1)-(N) (collectively referred to herein as 2650).

The architecture 2600 may further include an algorithm control block2660 that may be configured to generate various control signalsincluding the R/W (read/write) command to the memory 2610 and an MA(memory address) when needed. The number of available input signalsamples N_(s) may be divided into m_(s) segment each of length n_(s).During the first interval of n_(s)T_(s) seconds, the complex inputz_(c)(k) is connected to the complex to real converter that outputs adimension 2 vector z(k) with its elements equal to the real andimaginary parts respectively of the signal z_(c)(k). The signal z(k) isinput to the first stage AKF-MF 2630 that processes the n samples of thesignal z(k) to provide the estimates of various frequencies {circumflexover (ω)}_(m) (k) and frequency derivatives {circumflex over (α)}_(m)(k)for m=1, 2, . . . , N. The operation of the AKF-MF block 2630 has beenpreviously described with reference to FIG. 21.

Still referring to FIG. 26, the estimate of the m^(th) frequency{circumflex over (ω)}_(m)(k) and its derivative {circumflex over(α)}_(m)(k) may be provided to the correlation threshold detector CTDmblock 2640 for m=1, 2, . . . , N. The CTDm block 2640 may compute thecorrelation r_(m) ^(i) between the input signal and the complexexponential signal at frequency {circumflex over (ω)}_(m)(k) accordingto (104)-(108) and may compare the correlation with a predeterminedthreshold and further generate a signal t_(m) that is indicative of thecomparison at the output. The signal t_(m) may be equal to 1 if thethreshold condition is satisfied; otherwise it may be equal to 0. Thesignals t₁, t₂, . . . , t_(N) may, therefore, be provided to thealgorithm control block 2660. If any of the t_(m) values are equal to 0,the procedure is repeated with the stored n_(s) samples read from thememory 2610 and with the initial frequency estimates provided to thefirst stage AKF-MF 2630 selected in a random and independent manner froma selected distribution. For example, the initial frequency estimatesmay be selected from (−B/2 B/2) with a uniform distribution. For thefrequencies for which the threshold condition is satisfied, the initialfrequency estimates may be made equal to the corresponding estimatesobtained from the preceding run. This process of reading the n_(s)samples from the memory, processing them by the first stage AKF-MF 2630,evaluating correlations and comparing the correlation values with thethresholds in blocks CTD1, . . . , CTDN 2640 may be repeated until allt_(m) values are equal to 1 or the number of such trials exceeds somespecified value N_(t) which may, for example, be equal to 4. In analternative embodiment, instead of applying thresholds to individualcorrelation values, the threshold may be applied to the sum of thesquares of the individual correlation values evaluated in (110), and thethreshold may be some fraction of the total power P_(T) which may beeither known or estimated from the measurements as in (111). At thecompletion of the first process of the algorithm as described herein,the measurements may be processed by second stage AKFs 2650.

The samples of the input complex signal z_(c)(k) are inputted to a bankof stage 2 AKF1 to AKFN 2650 through the activation of switches S₁, S₂,and S₃ starting with the second segment of n_(s) samples in the firstrun of the second stage. Referring now to FIG. 27, a block diagram of anexample AKFm 2700 is illustrated and discussed. They AKFm 2700 may beused as the second stage AKFs 2650 of FIG. 26. The input complex signalz_(c)(k) may be multiplied by the output exp [jω_(m) ^(L)k] of anoscillator with frequency ω_(m) ^(L). The product of the complex inputand the oscillator output, z_(cm)(k) may be connected to a complex toreal converter 2720 that may output a dimension 2 vector z(k) with itselements equal to the real and imaginary parts respectively of thesignal z_(cm)(k). The signal z(k) is input to an AKF 2730 that processesthe n_(s) samples of the signal z(k) to provide estimates of variousfrequencies {circumflex over (ω)}_(m) ^(r)(k) and frequency derivatives{circumflex over (α)}_(m) ^(r)(k) for m=1, 2, . . . , N where{circumflex over (ω)}_(m) ^(r)(k) is the estimate of the differencebetween the frequency ω_(m)(k) and the oscillator frequency ω_(m) ^(L)that may be adjusted at the end of n_(s)T_(s) second intervals. Duringthe first run, i_(r)=1, of the second stage, the oscillator frequencymay be set equal to the estimate of the frequency ω_(m)(k) for k=n_(s),2n_(s), . . . , (m_(s)−1)n_(s) where m_(s) is the number of segments. Inthe second and subsequent runs of step 2, i_(r)>1, and for the case whenω_(m)(k) is constant over the period of N_(s) samples, the oscillatorfrequency at the beginning of the first segment is set equal to theestimate of ω_(m) obtained at the end of the previous run. The signalsI₁ to I₅, starting with initial values of zero, may change at the end ofn_(s)Ts second intervals by the clock signal CK2 (not shown) that may beincorporated in the algorithm control block 2660 of FIG. 26. In onenon-limiting example, I₁=1 if k=n_(s) and i_(r)=1, and ω_(m) ^(L) is setequal to the estimate of frequency {circumflex over (ω)}_(m) ⁰(k)obtained from the stage 1 AKF-MF at the end of processing by the stage 1AKF-MF as the output of the offset due to nonzero α_(m) is zero. As aresult, the outputs of the delay blocks in accumulators 1 and 2 of FIG.27 are equal to 0, and thus ω_(m) ^(L)={circumflex over (ω)}_(m) ⁰. Thesignal I₂ is equal to 1 if either k=i_(s)n_(s), i_(s)>1 or i_(r)>1otherwise I₂=0.

When I₂=1, the frequency estimate {circumflex over (ω)}_(m) ^(r)(k) maybe added on to the previous value of ω_(m) ^(L) providing an updatedvalue for ω_(m) ^(L). Thus ω_(m) ^(L) may be set to {circumflex over(ω)}_(m) ⁰(k) in the first run at k=n_(s) and from there on it may bemodified by {circumflex over (ω)}_(m) ^(r)(k) at intervals of n_(s)T_(s)seconds. The sum of {circumflex over (ω)}_(m) ^(r)(k) and ω_(m) ^(L) atthe end of a previous n_(s) length segment provides the frequencyestimates during any run and any segment. These estimates are timemultiplexed with the estimates {circumflex over (ω)}_(m) ⁰(k) from stage1 in MUX1 controlled by the algorithm control block details, which arenot shown. Thus if n_(r)=1, then the first segment final estimates areobtained from {circumflex over (ω)}_(m) ⁰(k) with the final estimatesfor the subsequent segments provided by the second stage. On the otherhand, if n_(r)>1, then all of the final estimates come from stage 2.

I₄ may be equal to 1 if k=n_(s) and i_(r)>1. Thus the signal ξ₁ may beequal to {circumflex over (α)}_(m) ⁰(k) if k=n_(s) and i_(r)=1; equal to{circumflex over (α)}_(m) ^(r)(k) if k=n_(s) and i_(r)>1; and equal to 0otherwise. An accumulator 2 2740 output may be set equal to ξ₁ atk=n_(s) for any value of i_(r). The signal I₅ is equal to 1 ifk=i_(n)n_(s), i_(s)>1. If I₅=1, the value {circumflex over (α)}_(m)^(r)(k) is added to the accumulated output. The signal I₃ is equal to 1if k=N_(s) and is 0 otherwise. At the end of N_(s) samples, anaccumulator 2740 output ξ₂ is the sum of the derivatives during them_(s) segments wherein it is assumed that during each segment thederivative is nearly constant. The accumulator 2750 output ξ₂ timesn_(s) is subtracted from ξ₃ to provide the oscillator frequency suchthat the oscillator frequency is nearly equal to the input signalfrequency at the beginning of n_(s) sample segments. The initialfrequency estimate for the AKF 2730 is reset to 0 at the end of eachn_(s) sample segment in each run, and the covariance matrix P(0/0) isreset to some specified value at the start of each run.

FIG. 28 illustrates an example block diagram of the correlationthreshold detector CTD m 2800. The CTD m 2800 may be used as the CTD m2640 of FIG. 26. The time-varying frequency estimates {circumflex over(ω)}_(m) ⁰(k) obtained during the last trial of the first stage may beinput to an accumulator 2810. The output of the accumulator {circumflexover (θ)}_(m)(k) representing the phase estimate may be input to an exp() block 2820 that may generate the complex exponential signals_(m)(k)=exp[−j{circumflex over (θ)}_(m)(k)] at the output. The signals(k) may be multiplied by the input complex signal z_(c)(k) in amultiplier 2830. The output of the multiplier 2830 may be averaged overan interval of n_(c) samples where n_(c) is a subinteger multiple ofn_(s) with n_(s)=m_(c)n_(c) and with m_(c) equal to an integer that maybe equal to 1. The output of an averaging block 2840 may be equal to thecoherent correlation r_(m) ^(c)(l) over the l^(th) interval of n_(c)samples according to (106). The coherent correlation may be provided toa cascade of an absolute value squared block 2850, an averaging block2860, and a square root computation block 2870, providing thesemi-coherent correlation r_(m) ^(i) over the interval of n_(s) samplesaccording to (107). The semi-coherent correlation r_(m) ^(i) may beinput to a threshold detector block 2880 that may compare the inputagainst a threshold value γ_(m) ^(th). If the input to the thresholddetector block 2880 exceeds the threshold r_(m) ^(th), then the outputt_(m) of the threshold detector block 2880 may be equal to 1. Otherwise,the output t_(m) of the threshold detector block 2880 may be equal to 0.The output t_(m) may be provided to the algorithm control block 2660 ofFIG. 26, which may determine the actions to be performed by the AFK-CDarchitecture 2600 based on the value of t_(m).

Simulation Results

FIGS. 29-33 show example plots of simulation results for the performanceof the AKF-CD algorithm for the case of a single frequency with thesignal-to-noise ratio of 10 dB. The total number of samples for thisnon-limiting example is 10. Both n_(s) and n_(c) are equal to 1. In thefirst phase, the value of N_(t) denoting the number of maximum trials is6. The frequency is selected randomly with a uniform probability densityfunction (pdf) in the interval (−200, 200) Hz with the sampling rateF_(s) equal to 400 Hz.

The amplitude A is equal to 1 volt. The correlation threshold γ_(m)^(th) in the first phase of the algorithm is 0.5. In the first trial ofthe first phase of the algorithm, the initial frequency estimate isselected to be 0 Hz; in the subsequent trials, the initial estimate ismade randomly with a uniform pdf in the range (−40, 40 Hz). The initialcovariance matrices are selected as P(0/0)=2 I₂ and P_(p)(0/0)=2 I₂ withI₂ equal to the 2×2 identity matrix. At the end of the first phase ofthe algorithm, the P_(p)(0/0) is reset to P_(p)(0/0)=0.05 I₂ to speed upthe convergence rate of the AKF. After that, it is updated according toAKF update equation (26).

The simulations were performed with 50 random frequency selections fromthe range of (−200, 200) Hz according to the uniform probability densityfunction (pdf). For each such frequency selection, the algorithm wassimulated with 100 independent noise sequences. FIG. 29 shows the RMSfrequency error in Hz obtained during the last trial of the first phaseof the algorithm for each of the 50 frequency selections. FIG. 30 plotsthe result when the AKF algorithm is repeated a second time afterfrequency shift according to (112). As may be inferred from thesefigures, for each frequency selection, the AKF-CD algorithm converges toa relatively small frequency estimation error. FIGS. 31 and 32 plot theRMS error obtained by averaging the errors obtained in the 50 frequencyselections. The results are also compared with the C-R bound given by(49a) for the same signal-to-noise ratio and the maximum frequency B=200Hz. As may be seen from FIG. 33, the difference between the AKF-CDalgorithm and the C-R bound is small. It may be emphasized that theinitial frequency estimates in FIGS. 31 and 32 are based on the completedata of 10 samples and, thus, it is the final error at the end of 10samples that is more meaningful in terms of the comparison with the C-Rbound. FIG. 33 shows the 50 frequencies that were selected in thesimulations. FIG. 34 shows the number of maximum trials in the firstphase of the algorithm for each of the 50 selected frequencies showingthat only 6 out of the 50 frequency selections required retrials forsome noise sequences.

FIGS. 35 and 36 plot the RMS frequency error averaged over all 50initial frequency selections for the case of SNR=7 dB and thefrequencies selected from the range of (−100, 100) Hz with N_(s)=20samples. All of the algorithm parameters are the same as those for the10 dB SNR case with n_(s)=n_(c)=10 and γ_(m) ^(th)=0.4.

FIG. 37 shows the RMS frequency error in Hz obtained during the secondrun of the algorithm for each of the 50 frequency selections for thecase of SNR=5 dB and the frequencies selected from the range of (−100,100) Hz with N_(s)=20 samples. All of the algorithm parameters are thesame as those for the 7 dB SNR case for the first run. The RMS error formost of the 50 cases lies between 2 Hz and 8 Hz. FIG. 38 plots the RMSfrequency error averaged over all 50 initial frequency selections forthe case of SNR=5 dB. The result for the second run (not shown) issimilar. In this case, although the error in absolute terms isrelatively small, there is a considerable difference between theestimation error obtained by the AKF-CD algorithm and the C-R bound.Further embodiments of the disclosure teach the AKF-FFT-CD algorithmwhich provides estimation errors close to the C-R bound for an SNR evensmaller than 0 dB.

Frequency Estimation for Very Low Signal-to-Noise Ratio Conditions

For relatively very low signal-to-noise ratio conditions, the Kalmanfilter based frequency estimation algorithm may not converge to thecorrect frequency. Instead it may converge to a frequency with anestimation error in the range of −F_(s)/2 to F_(s)/2 with some non-zeroprobability, i.e., the algorithm may diverge. Thus the algorithmconverges to the correct frequency for some noise sequence and theinitial estimate, while for some other noise sequence and/or initialfrequency estimates, it may converge to an estimate with a relativelylarge error with some non-zero probability. Moreover, the smaller thesignal-to-noise ratio, the relatively higher is the probability ofdivergence. For this case, another embodiment of the disclosure ispresented that provides good convergence even under very lowsignal-to-noise ratio conditions.

FFT-AKF-CD Algorithm

In this embodiment of the disclosure, the range of frequency uncertaintyF_(s)=2B is divided into M intervals of length Δf=F_(s)/M. The integer Mis selected such that the ratio Mγ=MA²/(2σ_(n) ²), with γ=A²(2σ_(n) ²)denoting the SNR, is higher than some specified threshold value γ_(th)and (N_(s)/M) is an integer. For example, with γ=0.25 and γ_(th)=5,M=20, let z_(c)(k) denote the complex valued signal asz _(c)(k)=z _(i)(k)+jz _(q)(k);j≡√{square root over (−1)}  (113)

The signal z_(c)(k) is correlated with complex exponential functions offrequency f_(m)=−B+mΔf; m=0, 1, . . . , M−1 over an interval ofM=F_(s)/Δf samples providing the correlation values

$\begin{matrix}{{{{r_{m}(0)} = {\frac{1}{M}{\sum\limits_{i = 1}^{M}\;{{z_{c}(i)}{\exp\left( {{- j}\; 2\pi\; f_{m}i\; T_{s}} \right)}}}}};{f_{m} = {{- B} + {m\;\Delta\; f}}};{m = 0}},1,\ldots\mspace{14mu},{M - 1}} & (114)\end{matrix}$

The expression in (114) can also be evaluated in terms of the discreteFourier transform (DFT) or the fast Fourier transform (FFT) bysubstituting the expression for f_(m) in the expression for r_(m)(0).Thus

$\begin{matrix}\begin{matrix}{{r_{m}(0)} = {\frac{1}{M}{\sum\limits_{i = 1}^{M}\;{{z_{c}(i)}{\exp\left( {j\; 2\;\pi\; B\; i\; T_{s}} \right)}{\exp\left( {{- j}\; 2\pi\; i\; m\;\Delta\;{fT}_{s}} \right)}}}}} \\{{{= {\frac{1}{M}{\sum\limits_{i = 1}^{M}\;{{z_{c}(i)}\left( {- 1} \right)^{i}\mspace{14mu}{\exp\left( {{- j}\; 2\pi\; m\text{/}M} \right)}{\exp\left\lbrack {{- \; j}\; 2{\pi\left( {i - 1} \right)}m\text{/}M} \right\rbrack}}}}};{m = 0}},1,\ldots,{M - 1}}\end{matrix} & (115)\end{matrix}$

From (115), the correlation value r_(m)(0) may be expressed as

$\begin{matrix}{{r_{m}(0)} = {\frac{1}{M}{\exp\left( {{- j}\; 2\pi\; m\text{/}M} \right)}{z_{0\; f}(m)}}} & (116)\end{matrix}$

where z_(0f)(f) is the discrete Fourier transform of z_(m)(i) given by

$\begin{matrix}{{{{z_{0\; f}(m)} = {\sum\limits_{i = 0}^{M - 1}\;{{z_{m}(i)}{\exp\left( {{- j}\; 2\pi\; i\; m\text{/}M} \right)}}}};{{z_{m}(i)} = {\left( {- 1} \right)^{i + 1}{z\left( {i + 1} \right)}}};{m = 0}},1,\ldots\mspace{14mu},{M - 1}} & (117)\end{matrix}$

In (117), the sequence z_(m)(i) is obtained by multiplying thealternative elements of the sequence z_(c)(i) by −1. The DFT of z_(m)(i)in (117) can also be evaluated by the M-point fast Fourier transform ofz_(m)(i).

With N₁=N_(s)/M, the correlation values are evaluated over thesubsequent (N_(I)−1) intervals by replacing the subsequence [z(1) z(2),. . . , z(M)] by the subsequence of z(k) for the respective interval.The estimate of the average power present during the interval of Nsamples of z_(c)(k) is obtained by averaging the absolute value squared|r_(m)|² over the N_(I) intervals for m=0, 1, . . . , M−1, i.e.,

$\begin{matrix}{P_{m} = \left. {\frac{1}{N_{I}}\sum\limits_{i = 0}^{N_{I} - 1}}\; \middle| {r_{m}(i)} \right|^{2}} & (118)\end{matrix}$

Substitution of z_(c)(i) from (3), (4), and (113) in (115) yields

$\begin{matrix}{{{r_{m}(0)} = {{s_{m}(0)} + {n_{m}(0)}}}{where}} & (119) \\{{{{s_{m}(0)} = {\frac{1}{M}{\sum\limits_{i = 1}^{M}\;{A\;{\mathbb{e}}^{j\;\varphi}{\exp\left( {j\; 2\pi\; f_{dm}i\; T_{s}} \right)}}}}};{f_{dm} = {f_{0} - f_{m}}};{m = 0}},1,\ldots\mspace{11mu},{M - 1}} & (120)\end{matrix}$

Evaluation of the sum in (120) results in

$\begin{matrix}{{{{s_{m}(0)} = {\frac{A\;{\mathbb{e}}^{j\varphi}}{M}{{\exp\left\lbrack {{j\left( {M + 1} \right)}\pi\; f_{dm}T_{s}} \right\rbrack}\left\lbrack \frac{\sin\left\lbrack {M\;\pi\; f_{dm}T_{s}} \right\rbrack}{\sin\left( {\pi\; f_{dm}T_{s}} \right)} \right\rbrack}}};{m = 0}},1,\ldots\mspace{14mu},{M - 1}} & (121)\end{matrix}$

For −Δf≦f_(dm)<Δf, Δf=1/(MTs) and M>>1, the absolute value of theexpression in (121) may be approximated by

$\begin{matrix}{{\left| {s_{m}(0)} \middle| {\cong A} \middle| \frac{\sin\left( {\pi\; f_{dm}\text{/}\Delta\; f} \right)}{\left( {\pi\; f_{dm}\text{/}\Delta\; f} \right)} \right|;{{{- \Delta}\; f} \leq f_{dm} < {\Delta\; f}};{m = 0}},1,\ldots\mspace{14mu},{M - 1}} & (122)\end{matrix}$

For the case whenf _(dm)=2nΔf+f _(δm) ;−Δf≦f _(δm) <Δf;n an integer with 0<|n|<M/4  (123)

The absolute value of s_(m)(0) in (122) may be approximated as

$\begin{matrix}\begin{matrix}\left| {s_{m}(0)} \middle| {\cong A} \middle| \frac{\sin\left( {\pi\; f_{\delta\; m}\text{/}\Delta\; f} \right)}{M\mspace{14mu}{\sin\left( {2\pi\; n\text{/}M} \right)}} \right| \\{{\left. {\cong \frac{\kappa\; A}{2\pi\; n}} \middle| {\sin\left( {\pi\; f_{\delta\; m}\text{/}\Delta\; f} \right)} \right|;{{{- \Delta}\; f} \leq f_{\delta\; m} < {\Delta\; f}};\left. {0 <} \middle| n \middle| {< \frac{M}{4}} \right.;{m = 0}},1,\ldots\mspace{14mu},{M - 1}}\end{matrix} & (124)\end{matrix}$

where κ is a constant in the range of 1 and 2. The noise term n_(m)(0)in (119) is given by

$\begin{matrix}{{{{n_{m}(0)} = {\frac{1}{M}{\sum\limits_{i = 1}^{M}\;{{n(i)}{\exp\left( {{- j}\; 2\pi\; f_{m}i\; T_{s}} \right)}}}}};{{n(i)} = {{n_{i}(i)} + {j\;{n_{q}(i)}}}};{j = \sqrt{- 1}};{m = 0}},1,\ldots\mspace{14mu},{M - 1}} & (125)\end{matrix}$

From (125) one obtains

$\begin{matrix}{{{\left| {n_{m}(0)} \right|^{2} = {{{n_{m}(o)}{n_{m}^{*}(0)}} = \left. {\frac{1}{M^{2}}\sum\limits_{i = 1}^{M}}\; \middle| {n(i)} \middle| {}_{2}{{+ \frac{2}{M^{2}}}{\sum\limits_{i = 1}^{M}\;{\sum\limits_{\substack{{k = 1}, \\ k > i}}^{M}\;{{Re}\left\{ {{n(i)}{n^{*}(k)}} \right\}}}}} \right.}};}{{{{n(i)} = {{n_{i}(i)} + {j\;{n_{q}(i)}}}};{j = \sqrt{- 1}};{m = 0}},1,\ldots\mspace{11mu},{M - 1}}} & (126)\end{matrix}$

In (126), Re( ) denotes the operation of taking the real part of theargument. As the sequence n(i) is zero mean and independent, it followsfrom (126) that

$\begin{matrix}{{{{E\left\{ \left| {n_{m}(0)} \right|^{2} \right\}} = \frac{2\sigma_{n}^{2}}{M^{2}}};{m = 0}},1,\ldots\mspace{14mu},{M - 1}} & (127)\end{matrix}$

wherein equation (127) E denotes the expected value operator. From (125)it follows that n_(m)(0) and n_(l)(0) are uncorrelated and thus due totheir Gaussian distribution are independent random variables as

$\begin{matrix}{{{{E\left\lbrack {{n_{m}(0)}{n_{k}^{*}(0)}} \right\rbrack} = {\frac{1}{M^{2}}{\sum\limits_{i = 1}^{M}\;{{E\left\lbrack \left| {n(i)} \right|^{2} \right\rbrack}{\exp\left\lbrack {j\; 2{\pi\left( {f_{k} - f_{m}} \right)}{\mathbb{i}}\; T_{s}} \right\rbrack}}}}};m},{k = 0},1,\ldots\mspace{14mu},{M - 1}} & (128)\end{matrix}$

For k≠m, the summation in (128) may be evaluated as

$\begin{matrix}\begin{matrix}{{{E\left\lbrack {{n_{m}(0)}{n_{k}^{*}(0)}} \right\rbrack} = {\frac{2\sigma_{n}^{2}}{M^{2}}{\sum\limits_{i = 1}^{M}\;{\exp\left\lbrack {j\; 2\pi\; f_{d}{iT}_{s}} \right\rbrack}}}};{f_{d} = {{f_{k} - f_{m}} = {\left( {k - m} \right)\Delta\; f}}}} \\{= {\frac{2\sigma_{n}^{2}}{M^{2}}{\exp\left\lbrack {{j\left( {M + 1} \right)}\pi\; f_{d}T_{s}} \right\rbrack}\frac{\sin\left( {\pi\;{Mf}_{d}T_{s}} \right)}{\sin\left( {\pi\; f_{d}T_{s}} \right)}}} \\{= {\frac{2\sigma_{n}^{2}}{M^{2}}{\exp\left\lbrack {{j\left( {M + 1} \right)}\pi\; f_{d}T_{s}} \right\rbrack}\frac{\sin\left\lbrack {\pi\left( {k - m} \right)} \right.}{\sin\left\lbrack {{\pi\left( {k - m} \right)}\text{/}M} \right\rbrack}}}\end{matrix} & (129)\end{matrix}$

The last term in (129) is zero as |k−m|≦M−1 resulting inE[n _(m)(0)n* _(k)(0)]=0;k≠m  (130)

Taking the absolute value squared on both sides of (119) yields|r _(m)(0)|² =|s _(m)(0)|² +|n _(m)(0)|²+2Re{s* _(m) n _(m) };m=0,1, . .. ,M−1  (131)

From (125) to (130), the expression in (131) may be expressed as

$\begin{matrix}{{{\left| {r_{m}(0)} \right|^{2} = \left| {s_{m}(0)} \middle| {}_{2}{{+ \frac{2\sigma_{n}^{2}}{M}} + {\xi_{m}(0)}} \right.};{m = 0}},1,\ldots\mspace{14mu},{M - 1}} & (132)\end{matrix}$

In (132) ξ_(m)(0) for m=0, 1, . . . , M−1 are zero mean independentrandom variables, and s_(m)(0) is given by (122). The average power atfrequency f_(m) is estimated by

$\begin{matrix}{{{P_{m} = {\left. {\frac{1}{N_{I}}\sum\limits_{i = 0}^{N_{I} - 1`}}\; \middle| {r_{m}(i)} \right|^{2} = \left| {s_{m}(0)} \middle| {}_{2}{{+ \frac{2\sigma_{n}^{2}}{M}} + \xi_{a,m}} \right.}};{\xi_{a,m} = {\frac{1}{N_{I}}{\sum\limits_{i = 0}^{N_{I} - 1}\;{\xi_{m}(i)}}}};{m = 0}},1,\ldots\mspace{14mu},{M - 1}} & (133)\end{matrix}$

In (133), ξ_(a,m) for m=0, 1, . . . , M−1 are zero mean independentrandom variables with their variance decreasing inversely withincreasing N_(I). The unknown frequency f₀ may be estimated bymaximization of P_(m) over m. However, for the case of low SNR γ, andlimited number of samples N and thus small N_(I), ξ_(a,m) hassignificant variance resulting in a significant probability of selectingan incorrect value of m and thus causing a relatively large frequencyestimation error f₀−f_(m). Moreover, even when the value of m iscorrectly selected, there is an estimation error in the range of −Δf/2and −Δf/2. From (133) the value of |s_(m) (0)|² may be onlyA²sinc²(0.5)=0.4 A² when f₀ is in the middle of any two of thefrequencies f_(m) resulting in about a 4 dB loss in signal powercompared to the case when the frequency f₀ is equal to one of thefrequencies f_(m) where f_(m) is given by (114). Therefore, theprobability of selecting an incorrect frequency f_(m) is much higher inthe case when the unknown frequency falls in between the two frequenciesf_(m) compared to the latter case.

In order to minimize the probability of selecting an incorrect frequencyf_(m), a number N_(f) of frequencies {f₀ ^(s), f₁ ^(s), . . . , f_(N)_(f) ₋₁ ^(s)} are selected from the set of frequencies {f₀, f₁, . . . ,f_(M-1)} defined in (114) with the highest average power P_(m) definedin equation (133). The adaptive Kalman filter algorithm described inequations (12)-(36) is repeated N_(f) times with the stored in-phase andquadrature observations z_(i)(k) and z_(q)(k) defined in equations (3)and (4) with the initial frequency estimates equal to the f_(m) ^(s),m=0, 1, . . . , N_(f)−1 and with reduced initial prediction errorcovariance matrix P_(p)(0/0) reduced by a factor of the order M². Theinitial covariance matrix is selected equal to εI₂ with I₂ denoting the2×2 identity matrix. Thus ε is the variance of the error in theestimates cos [2π(f0+δf)Ts] and sin [2π(f0+δf)Ts] for the parameters β₁and β₂ respectively, where f₀ denotes the true frequency and δf denotesthe error in the frequency estimate with its variance equal to σ_(f) ².With the approximationcos [2π(f ₀ +δf)Ts]−cos (2πf ₀ Ts)≅−2πT _(s) δf sin (2πf ₀ Ts)  (134)

the estimate for ε is given by

$\begin{matrix}{ɛ = {{\frac{4\pi^{2}}{F_{s}^{2}}\sigma_{f}^{2}{E\left\lbrack {\sin^{2}\left( {2\pi\; f_{0}T_{s}} \right)} \right\rbrack}} = {\frac{2\pi^{2}}{F_{s}^{2}}\sigma_{f}^{2}}}} & (135)\end{matrix}$

The second equality in (76) follows from the fact that the frequency f₀itself may have any value in the interval (−F_(s)/2, F_(s)/2) and thusthe argument of the sin function may lie in the interval (−π, π)resulting in the average for sin²( ) equal to 0.5. If the frequencyerror δf can have any value in the interval (−Δf/2, Δf/2), then σ_(f) ²may be approximated by ⅓(Δf/2)²=(Δf)²/12. Substitution for σ_(f) ² in(135) results in the estimate for ε equal to π²/(6M²) or approximately1/m². For M=20, this estimate is equal to 0.0025. This is much smallercompared to the value for P_(p)(0/0)=4 I₂ used in the EKF without theFFT step. The estimate of the ε derived in (135) applies to thatparticular initial frequency estimate f_(m) ^(s) that is close to thetrue frequency f₀ with the error smaller than Δf/2 in its magnitude. Theestimate for P_(p)(0/0) is based on the initial knowledge and needs tobe only approximate for the AKF to converge. Denoting by f_(m) ^(K),m=0, 1, . . . , N_(f)−1 the final frequency estimates obtained by theAKF corresponding to the N_(f) initial frequency estimates f_(m) ^(s),m=0, 1, 2, . . . , N_(f)−1, an estimate of the power P_(m) ^(K) isobtained similar to that for the P_(m) in equation (133). Thus,

$\begin{matrix}{{{P_{m}^{K} = \left. {\frac{1}{N_{I}}\sum\limits_{i = 0}^{N_{I} - 1}}\; \middle| {r_{m}^{K}(i)} \right|^{2}};{m = 0}},1,\ldots,{M - 1}} & \left( {136a} \right) \\{{{{r_{m}^{K}(i)} = {\frac{1}{M}{\sum\limits_{k = 1}^{M}\;{{z_{c}\left\lbrack {k + {\left( {i - 1} \right)M}} \right\rbrack}{\exp\left( {{- j}\; 2\pi\;{kf}_{m}^{K}T_{s}} \right)}}}}};{j = \sqrt{- 1}};{m = 0}},1,\ldots,{{N_{f} - 1};{i = 0}},1,\ldots,{N_{I} - 1}} & \left( {136b} \right)\end{matrix}$

In the absence of any initial estimate of the received signal powerlevel A², the frequency f_(m) ^(K) corresponding to the highest value ofP_(m) ^(K) is selected for further processing. However, with someestimate of the signal power level A², detection is performed on thehighest value of P_(m) ^(K) to eliminate the possibility that none ofthe frequencies f_(m) ^(K) are close to the true frequency f₀. Theprobability of this is a function of M, N_(I), Nf, and the SNR γ and maybe indicated as the highest power out of the M values P_(m) ^(K) beingrelatively very small.

For the frequency f_(m) ^(K) close to f₀, the expected value of thecorresponding power P_(m) ^(K) obtained from (136) or from the equations(121) and (133) with f_(dm) in (121) equal to 0 will be close toP_(a)=A²+2σ_(n) ²/M. The actual power level may have mean P_(a) with anerror with a variance that is a function of N₁ and γ. The actual powerwill exceed the value P_(a) with a probability of 0.5. Thus a thresholdvalue P_(T) on the power level P_(m) ^(K) can determine whether or notthe corresponding frequency f_(m) ^(K) is close to f_(m). However, ifP_(T) is too high, one may miss detecting the correct frequency f_(m)^(K) whereas lowering the threshold to a relatively small value mayresult in a false detection. As an example, with A²=1, σ_(n) ²=2.5corresponding to the SNR γ=0.2, and with M=20, one obtains P_(a)=1.25and a threshold value of P_(T)=0.8 may be selected. If the highest valueof P_(m) ^(K), m=0, 1, . . . , M−1, does not exceed P_(T), then the sumof the two highest powers P_(m) ^(K) is compared with the thresholdP_(T2). If the sum does exceed the threshold P_(T2) and thecorresponding two frequencies are within Δf/2 of each other, then thefrequency with the highest power P_(m) ^(K) among the two frequencies isaccepted as the correct intermediate estimate of the frequency f₀.Alternatively the average value of the two frequencies may be taken asthe correct intermediate estimate of the frequency f₀. For the examplewith P_(a)=1.25, a value for P_(T2) equal to 1 is selected. However, ifnone of the two threshold conditions are satisfied then the receivedsequence z_(c)(k) is further processed as follows.

For a specified number of trials N_(r)=2^(κ) for some integer κ≧0, thesteps of correlating the received sampled signal z_(c)(k) with thecomplex exponential signals with M frequencies f_(m)=−B+mΔf; m=0, 1, . .. , M−1, selecting N_(f) of these frequencies with the highestcorrelation values, running the AKF with each of these N_(f) frequenciesas the initial frequency estimates, correlating z(k) with the N_(f)final frequency estimates provided by the AKF, and comparing thecorrelation values with the two thresholds, are repeated until thethreshold on the correlation value is satisfied or the number of suchruns exceeds N_(r-1). In each such run, the M frequencies f_(m) definedin (114) are offset by a different value of δf. The offset values of δfin different runs are selected from the set of values {i Δf/N_(r); i=1,2, . . . , N_(r)−1}. However, the sequence in which these offset valuesmay be taken may not be in an increasing order. For example, withN=2^(κ) for an integer κ, the values may be taken in the order [2^(κ-1),2^(κ-2)×(1, 3), . . . , (1, 3, . . . , N_(r)−1)]Δf/N_(r). For examplewith N_(r)=8, the order is [4, 2, 6, 1, 3, 5, 7]Δf/8 or

$\left\lbrack {\frac{1}{2},\frac{1}{4},\frac{3}{4},\frac{1}{8},\frac{3}{8},\frac{5}{8},\frac{7}{8}} \right\rbrack\Delta\;{f.}$

The correlation P_(m) ^(d) between z_(c)(k) and the complex exponentialsignals with M frequencies f_(m) ^(d)=−B+mΔf+δf; m=0, 1, . . . , M−1,can be obtained in terms of generalized FFT as for the case of δf=0.

$\begin{matrix}{{{P_{m}^{d} = \left. {\frac{1}{N_{I}}\sum\limits_{i = 0}^{N_{I} - 1}}\; \middle| {r_{m}^{d}(i)} \right|^{2}};{m = 0}},1,\ldots\mspace{14mu},{M - 1}} & (137) \\{{r_{m}^{d}(i)} = {\frac{1}{M}{\exp\left( {{- j}\; 2\pi\; m\text{/}M} \right)}{z_{if}^{d}(m)}}} & (138) \\{{{{{z_{if}^{d}(m)} = {\sum\limits_{k = 1}^{M}\;{{z_{c}^{d}\left\lbrack {k + {\left( {i - 1} \right)M}} \right\rbrack}{\exp\left( {{- j}\; 2\;\pi\;{km}\text{/}M} \right)}}}};{i = 0}},1,\ldots\mspace{14mu},{N_{I} - 1}}{{{{z_{c}^{d}(k)} = {\left( {- 1} \right)^{k + 1}{z_{c}\left( {k + 1} \right)}{\exp\left\lbrack {{- \; j}\; 2\pi\;{kl}\text{/}\left( {MN}_{r} \right)} \right\rbrack}}};{{\delta\; f} = \left( {l\;\Delta\; f\text{/}N_{r}} \right)};{m = 0}},1,\ldots\mspace{11mu},{M - 1}}} & (139)\end{matrix}$

In equation (139), z_(if) ^(d)(m) m=0, 1, . . . , M−1, is the DFT of thei^(th) segment of length M of the modified signal z_(c) ^(d)(k) obtainedfrom z_(c)(k) as shown in equation (139) and thus can be evaluated by anFFT algorithm.

Denoting by f_(I) the estimate obtained by the combined FFT and EKFalgorithms, a further refinement in the estimate is achieved bydemodulating the signal z_(c)(k) by the complex exponential signal atfrequency f_(I) which is equivalent to changing the frequency of thesignal to f_(r)=(f₀−f_(I)), averaging the resulting signal over segmentsof length M, and processing the averaged signal by the AKF. Thus

$\begin{matrix}{{{{z_{s}(k)} = {{z_{c}(k)}{\exp\left\lbrack {{- \; j}\; 2\pi\;{kf}_{I}\text{/}F_{s}} \right\rbrack}}};{k = 1}},2,\ldots\mspace{14mu},N_{s}} & (140) \\{{{{z_{a}(k)} = {\frac{1}{M}{\sum\limits_{i = 1}^{M}\;{z_{s}\left\lbrack {i + {\left( {k - 1} \right)M}} \right\rbrack}}}};{k = 1}},2,\ldots\mspace{14mu},{N_{I} - 1}} & (141)\end{matrix}$

In (141), the signal z_(s)(k) is averaged over M consecutive samplesthus reducing the sampling rate to F_(sa)=F_(s)/M. This is in view ofthe fact that in the absence of a detection error, the frequencyestimation error in f_(I) is limited in magnitude to Δf/2 instead ofF_(s)/2. Thus the sampling rate is reduced by a factor of F_(s)/Δf=M.The averaging process reduces the noise variance also by the factor M.However, the normalized frequency uncertainty [(f₀−f_(I))F_(sa)] lies inthe range of (−π, π) and therefore the initial covariance matrix P(0/0)may be set equal to ε_(a)I₂ with ε_(a)=(π/3)²≅1 assuming a Gaussiandistribution for the initial frequency estimation error with 3σ valueequal to π. In the AKF, the sampling time T_(s) is replaced byT_(sa)=MT_(s). The final estimate of the unknown frequency f₀ is givenby{circumflex over (f)} ₀ =f _(I) +{circumflex over (f)} _(r)  (142)

where {circumflex over (f)}_(r) denotes the estimate of the frequencyf_(r)=(f₀−f_(I)) by the AKF. Finally accurate estimates of the amplitudeA and phase φ are obtained by correlating z_(c)(k) with the complexexponential function of frequency {circumflex over (f)}₀. Thus theestimates of the amplitude A denoted by Â and that of φ denoted by{circumflex over (φ)}, respectively, are given by

$\begin{matrix}{{{\hat{A} = |r|};{\hat{\varphi} = {\arg(r)}}}{where}} & (143) \\{r = {\frac{1}{N_{s}}{\sum\limits_{k = 1}^{N_{s}}\;{{z_{c}(k)}{\exp\left\lbrack {{- \; j}\; 2\pi\; k{\hat{f}}_{0}T_{s}} \right\rbrack}}}}} & (144)\end{matrix}$

FIG. 39 is a block diagram illustrating an architecture 3900 accordingto another embodiment of the disclosure for the case of very lowsignal-to-noise ratio. The input complex signal may be connected to theinput of a memory block 3910 and to a switch S₁. The output of theswitch can be selected either directly from the input or read from thememory 3910. The output of the switch S₁ may be connected to the inputof switch S₂ and depending upon the state of switch S₂ and S₃ it may berouted to stage 1 estimator 3920, stage 2 estimator 3930, or stage 3estimator 3940. With the position of switches S₁ and S₂ as depicted inFIG. 39, the input signal z_(c)(k) may be shifted in frequency byδf=1Δf/N_(r) by a frequency converter 3916 with N_(r)=2^(κ) for aninteger κ>0, 1=0, 1, . . . , N_(r)−1, Δf=F_(s)/M, and M is the FFT size.The output of the frequency converter 3916 is multiplied by (−1)^(k).The modified signal z_(m)(k) is input to serial to parallel S/Pconverter 3912 that splits the input signal z_(m)(k) into M outputs. TheM outputs of the S/P converter are input to the M-point FFT block 3918that evaluates the M-point FFT transform of the input and provides theFFT transform [z_(if)(0) z_(if)(1) . . . z_(if)(M−1)] at the output. Forthe case of i^(th) input block of size M, and for the case of 0frequency shift, δf=0, the operation on z_(c)(k) resulting in the FFTtransform is described by (117). The M outputs of the FFT block aremultiplied by (1/M) exp [−j2πm/M]; m=0, 1, . . . , M−1 in multipliers3919 and 3921 according to (116) illustrated for the case of i=0^(th)block, providing the correlation values r₀(i), . . . , r_(M-1)(i) at theoutputs of multipliers.

When the correlation values are input to a cascade of mod square andaveraging blocks 0, 1, . . . , M−1 3926 to 3936 providing the estimatedaverage power P₀, P₁, . . . , P_(M-1) present during the N_(s) samplescorresponding to N_(I) blocks of M samples each in the respectivepassbands of the M FFT filters, the FFT transform may be treated as abank of M FFT filters. Referring to (114)-(117), the m^(th) output ofthe FFT block 3918 corresponds to shifting the input signal frequency byf_(m)=−B+mΔf, and averaging the frequency shifted signal over aninterval of M samples. The two operations of frequency shifting andaveraging are equivalent to bandpass filtering the signal with a complexbandpass filter of center frequency f_(m) and zero crossing bandwidthequal to 2Δf, and down converting the filtered signal to baseband. Thepower levels P₀, P₁, . . . , P_(M-1) at the output of M averaging blocksare input to the frequency selector 1 3938 block that selects N_(f)frequencies from the set of frequencies f_(m), m=0, 1, . . . , M−1 withthe highest power levels. The selected frequencies {f₀ ^(s), f₁ ^(s), .. . , f_(N) _(f) ₋₁ ^(s)} are input to the respective AKF subsystemsAKFSS0, . . . , AKFSS(N_(f)−1) of the stage 2 estimator 3930 as theinitial frequency estimates for the respective AKFSS 0, . . . , N_(f)−1.

The block diagram of the adaptive Kalman filter subsystem (AKFSS) asdiscussed in reference to FIG. 2 may include the matrices given byQ(k−1)=ε_(Q)I₂ and R(k)=σ²I₂, σ²=σ_(n) ²/M, Q_(p)(k−1)=ε_(Qp)I₂ withε_(Qp) possibly equal to 0, and F_(p)(k−1)=I₂. With continuing referenceto FIG. 39, the input signal z_(c)(k) may be read from the memory 3910and provided to a complex to real converter 3914 to generate signal z(k)at the output. The signal z(k) may be input to the AKF subsystems AKFSS0, . . . , N_(f)−1 that also receive the frequencies {f₀ ^(s), f₁ ^(s),. . . , f_(N) _(f) ₋₁ ^(s)} as the initial frequency estimates with theinitial parameter error covariance matrix to all the AKFSSs given byP_(p)(0/0)=(1/M²) I₂. The input signal z(k), k=1, 2, . . . , N_(s) isprocessed by the (N_(f)−1) AKFSSs and may generate the final estimatesof frequencies f₀ ^(K), f₁ ^(K), . . . , f_(N) _(f) ₋₁ ^(K),respectively. The reset input I_(R1) from the algorithm control block3970 may provide the initialization for the AKFSS. At the end ofprocessing by the AKFSSs, the frequency estimates f₀ ^(K), f₁ ^(K), . .. , f_(N) _(f) ₋₁ ^(K) may be provided to correlator blocks corr 0, . .. , corr (N_(f)−1) blocks which are selected via the switch S₄ from thememory 3910. The correlator blocks corr 0, . . . , corr (N_(f)−1) mayestimate the average power present in the input signal at thefrequencies f₀ ^(K), f₁ ^(K), . . . f_(N) _(f) ₋₁ ^(K), respectively.

The frequencies f₀ ^(K), f₁ ^(K), . . . , f_(N) _(f) ₋₁ ^(K), along withthe estimates of the average power levels P₀ ^(K), P₁ ^(K), . . . ,P_(N) _(f) ₋₁ ^(K) present at these frequencies, are provided to thefrequency selector 2 3950. The frequency selector 2 3950 selects andoutputs the two frequencies f_(m) ₀ ^(K) and f_(m) ₁ ^(K), with highestaverage power levels P_(m) ₀ ^(K) and P_(m) ₁ ^(K), respectively, fromthe N_(f) input frequencies. The highest power level P_(m) ₀ ^(K) may beinput to the threshold detector 1 3954(1) with associated thresholdlevel P_(T1). The output of the threshold detector t₀ may be equal to 1if P_(m) ₀ ^(K)≧P_(T) ₁ and may be 0 otherwise. The sum of P_(m) ₀ ^(K)and P_(m) ₁ ^(K) may be input to the threshold detector 2 3954(2) withassociated threshold level P_(T2). The absolute value of the differencebetween f_(m) ₀ ^(K) and f_(m) ₁ ^(K) may be input to the thresholddetector 3 3954(3) where it may be compared to the resolution frequencyf_(res) and the output of the threshold detector may be 1 if |f_(m) ₀^(K)−f_(m) ₁ ^(K)|≦f_(res) and 0 otherwise. The outputs of thresholddetectors 2 and 3 3954(2), 3954(3) may be input to an AND gate 3960 thathas output t₁ equal to 1 if both the inputs are 1 and 0 otherwise. Boththe signals t₀ and t₁ may be provided to the algorithm control block3970. The frequency f_(m) ₀ ^(K) and the average frequency 0.5(f_(m) ₀^(K)+f_(m) ₁ ^(K)) may be input to a frequency selector 3 3976. On thecommand I₅ from the algorithm control block 3970, the frequency selector3 3976 may output one of the two frequencies at the output denoted byf^(I).

Based on the values of t₀ and t₁ the algorithm control block 3970 maygenerate the signal I₅ which has 3 possible values of 0, 1, and 2. Basedon I₅, the frequency selector 3 selects frequency f_(m) ₀ ^(K) if t₀=1,selects frequency 0.5(f_(m) ₀ ^(K)+f_(m) ₁ ^(K)) if t₀=0 and t₁=1, andselects none if both t₀ and t₁ are equal to 0.

If either t₀=1 or t₁=1, the first and second steps of the estimator arecomplete. However if both t₀ and t₁ are equal to 0, the algorithmcontrol block 3970 may select a different value of 1 or δf=1Δf/N_(r) andthe estimation procedure is repeated until either the thresholdcondition is satisfied or all of N_(r) values of δf have been processed.The frequency estimate f^(I) at the completion of step 1 and 2 may betaken as the final estimate of the unknown frequency f₀. However, afurther refinement in the frequency estimate may be made by the stage 3of the estimator block 3940.

FIG. 40 illustrates a block diagram of an example estimator 4000. Theestimator 4000 may be implemented as the stage 3 estimator 3940 of FIG.39. The complex signal z_(c)(k) is input to a frequency shifter thatshifts the frequency of the signal z_(c)(k) by f^(I). The frequencyestimate f^(I) at the output of the stage 2 of the estimator 3930 ofFIG. 39 may be inputted to an oscillator 4010 of FIG. 40 to set thefrequency of the oscillator to f^(I). The output of the oscillator 4010may be connected to a first frequency shifter 4014. The output z_(s)(k)of the first frequency shifter 4014 given by (140) may be input to anaveraging circuit 4020 that is configured to average the input z_(s)(k)over an interval of M samples and generates at its output an averagesignal z_(ac)(k) according to (141). The output of the averaging circuit4020 may be input to a complex to real converter 4026 that may beconfigured to convert the complex signal into a real vector signalz_(a)(k) that may be inputted to an AKF block 4030. The sample time tothe AKF block 4030 may be set to T_(sa)=MT_(s) with the initialestimates for the frequency ω⁰, and the parameter estimation errorcovariance matrix P(0/0) is set to 0 and I₂, respectively. Themeasurement error covariance matrix R(k) is set equal to (σ_(n) ²/M)I₂.The AKF block 4030 may process the N_(s)/M samples of the signalz_(a)(k) to generate the frequency estimate {circumflex over (f)}_(r)(k)that is the estimate of the error (f₀−f^(I)). Adding the frequencyestimate f^(I) from stage 2 of the estimator to {circumflex over(f)}_(r)(k) provides the final estimate {circumflex over (f)}₀ for theunknown frequency f₀.

{circumflex over (f)}₀ may control the frequency of a second oscillator4040 that generates the complex exponential signal at frequency{circumflex over (f)}₀ which may be input to the second frequencyshifter 4050. The other input of the second frequency shifter 4050 maybe with the complex signal z_(c)(k). The output of the second frequencyshifter 4050 may be inputted to an averaging circuit 4054 that mayaverage out the N_(s) samples as the output of the second frequencyshifter 4050 and generate the output r according to (145). The output ofthe averaging circuit 4054 may be inputted to the absolute value ∥ andthe arg( ) blocks, respectively, to provide the estimate of theamplitude Â and phase {circumflex over (φ)}, respectively of the inputsignal. When the expected RMS frequency estimation error is notinsignificant as may be assessed from the C-R bound, instead of usingthe coherent averaging in the estimation of Â, a quasi coherentcorrelation may be used in the estimation of the amplitude A. A switch(not shown) at the input of stage 3 estimator 4000 may route the inputz_(c)(k) to either the first or second frequency shifter 4014, 4050 asneeded by the algorithm control block 3970 of FIG. 39.

Referring now to FIG. 41 a block diagram of an example correlation block4100 according to embodiments of the disclosure is discussed. The corr mblocks of FIG. 39 may be implemented as correlation block 4100 of FIG.41. The complex signal z_(c)(k) may be inputted to a frequency shifter4110. The other input of the frequency shifter 4110 may be connected tothe output of an oscillator 4120. The oscillator 4120 frequency may beset equal to f_(m) ^(K) as provided by the AKSSS m of FIG. 39. Theoutput of the frequency shifter 4110 may be input to a first averagingcircuit 4130 which may be configured to average the input over aninterval of M samples generating r_(m) ^(K) at its output. The outputsignal r_(m) ^(K) may be inputted to a ∥² block 4140 that provides theabsolute value squared of the input. The ∥² block 4140 output|r_(m)^(K)|² may be inputted to a second averaging circuit 4150 that averagesout its input over the N_(I) samples providing the output P_(m) ^(K)according to (137a). The reset input I_(R2) from the algorithm controlblock 3970 of FIG. 39 may provide the initialization for the averagingcircuits 4130, 4150.

Estimation of Multiple Frequencies with the AKF-EFT-CD Algorithm

The AKF-EFT-CD algorithm can be applied to the case of the input signalthat includes N signals of different frequencies as in (74). FIG. 42shows the block diagram of the FFT-AKF-CD estimator 4200 for theestimation of multiple frequencies. In a first stage 4210 of theestimator, the interval Δf in (114) may be selected so that the minimumspacing among the input frequencies present in the input signal ishigher than 2Δf and thus M=F_(s)/Δf is higher than 2N. As for the caseof single frequency, the average power P_(m) may be estimated from (114)to (118) for m=0, 1, . . . , M−1. The number of frequencies N_(f) withthe highest power P_(m) may be selected for the AKF estimation withN≦N_(f)≦M, wherein a relatively higher value of N_(f) may result in arelatively small probability of missing any of the input frequencies inthe set of N_(f) selected frequencies. The set of frequencies thusselected are denoted by f_(m) ^(s), m=0, 1, . . . , N_(f)−1.

In a second stage 4220 of the estimator, the AKF algorithm is applied tothe signal model described in (74) to (87) with N replaced by N_(f)therein. In the signal model, the (N_(f)−N) of the amplitudes A_(m) arezero as the number of frequencies N_(f) selected is higher than N. Withf_(m) ^(K), m=0, 1, . . . , N_(f)−1 denoting the final frequencyestimates obtained by the AKF corresponding to the N_(f) initialfrequencies f_(m) ^(s), the average power P_(m) ^(K) at thesefrequencies is evaluated by equation (136). The N frequencies with thehighest power estimates P_(m) ^(K) are selected as the frequencyestimates are obtained. However, if any two frequency estimates arewithin the resolution frequency f_(res) of each other, with f_(res)equal to a predetermined fraction of Δf, then these two estimates areconsidered to be the estimates of a single frequency out of the Nfrequencies. The two such frequency estimates are averaged to providethe overall estimate of the frequency with the corresponding powerestimates added to obtain the overall estimate of the power at thatfrequency. After such a combining of the powers, the N frequencies withthe highest power estimates are selected as the N frequency estimatesform the step 2 of the algorithm.

As an option, a threshold on the estimates of the power levels of theindividual frequencies or a threshold on the total power in all of the Nfrequencies may be performed as for the case of single frequencyestimation. If the threshold condition is not satisfied, then theprocedure may be repeated with the M frequencies f_(m)=B+mΔf, m=0, 1, .. . , M−1 offset by 1Δf/2^(κ) for a selected integer κ. with the offsetinteger taking possible values 1, 2, . . . , 2^(κ)−1, as for the case ofthe single frequency estimation, to obtain a new set of initialfrequencies f_(m) ^(s), m=0, 1 . . . , N_(f)−1, running the AKFalgorithm with these new estimates and obtaining the frequency estimatesf_(m) ^(K), m=0, 1, . . . , N_(f)−1 and selecting N frequencies out ofthe N_(f) estimates based on their power estimates. The procedure may berepeated with different values of the offset integer i until thethreshold condition is satisfied or all the 2^(κ)−1 values of 1 areexhausted.

In the third stage 4230(1)-(N) estimator, the estimates f_(m) ^(I) ofeach of the N frequencies obtained from the second stage 4220 arefurther individually refined. In this step, the signal z_(c)(k) may bemultiplied by complex exponential function of frequency f_(m) ^(I) andthe product averaged over an interval of M=F_(s)/Δf samples as in(140)-(141) to obtain the averaged signal z_(a,m)(k) that is applied tothe AKF algorithm with the dimension of the state vector equal to 2 toobtain the estimate {circumflex over (f)}_(m) ^(r) of the residualfrequency f_(m) ^(r) that is equal to the error in the frequencyestimate. The final frequency estimate is given by (145) similar to theestimate in (142) for the single frequency case.{circumflex over (f)} _(m) =+f _(m) ^(I) +{circumflex over (f)} _(m)^(r) ,m=1,2, . . . N  (145)

In the estimation of multiple frequencies, the signal z_(c)(k) may besegmented into multiple signals by a polyphase filter or the bank of FFTfilters 4240 as in (94-96), and the AKF-FFT-CD algorithm is applied toeach of these segments to reduce the dimension of the state vector inthe AKF algorithm and thus the computational requirements.

Akf-Fft-Cd Algorithm Without the Last Step of Frequency Translation andDecimation

Simulation results are presented first for the case when the refinementstep of frequency translation and decimation in equations (140) and(141) followed by the application of the AKF to the resulting signal isnot performed. Instead the estimate based on selecting one of the N_(f)estimates obtained by the AKF on the basis of the correlation detectionperformed on these estimates is used as the final estimate. FIG. 43shows an example in-phase component of the received complex signalz_(c)(k) with amplitude 1 volt and frequency of 121.2 Hz with an SNR ofγ=−5 dB that appears like noise because the noise power is 3 timeshigher than the signal power. The AKF-FFT-CD algorithm is applied forthe estimation of the frequency and amplitude with the number ofsimulation runs equal to 100 with the initial covariance matrix equal to4I₂ as in the previous examples corresponding to the frequency range of[−200, 200] Hz. The value of M is selected equal to 20. FIG. 44 shows aplot of the frequency estimation error at the end of 80 samples for the100 simulation runs corresponding to a cumulative SNR defined asγ_(c)=10 log(N_(s)γ) of 14 dB. The RMS frequency error from the 100simulation runs is equal to 1.23 Hz compared to 0.4 Hz from the C-Rbound. FIG. 45 plots the amplitude estimates obtained from equation(143) for the 100 simulation runs.

FIG. 46 plots the result for the case of SNR=−7 dB after processing 160samples. The signal frequency selected randomly was −160.5 Hz for thiscase. The initial covariance matrix P_(p)(0/0) is equal to 2I₂ and theRMS frequency error is 0.70 Hz. However, if the two outliers are notincluded in the RMS error, the RMS error is only 0.32 Hz compared to theC-R bound of 0.17 Hz.

FIG. 47 plots the magnitude of the frequency estimation error for thecase of SNR=−10 dB after processing 200 samples corresponding to γ_(c)equal to 13.0 dB. The signal frequency selected randomly was −138.5 Hzfor this case. The RMS frequency estimation error for this case is 1.07Hz after rejecting 1 outlier at 325 Hz or equivalently at −75 Hz. TheRMS frequency error can be further reduced to bring it in line with theC-R bound by increasing the value of N_(f) and/or N_(r). However, inview of the relatively small value of the resulting RMS error, increasein the values of N_(f) or N_(r) may not be necessary for this case. Forstill lower SNRs, for example, at an SNR of −12 dB even after processing500 samples corresponding to γ_(c) equal to 15 dB, there may be 5outliers as shown in FIG. 48, the modulo F_(s) estimation error for theoutliers as lies in the range of about 10 Hz to 90 Hz. However,increasing N_(f) to 10 from 5 and N_(r) to 8 from 4 reduces the numberof outliers is reduced to 2 as shown in FIG. 49. Further increase inN_(f) and/or N_(r) may remove the remaining two outliers. FIG. 50 showsthe estimation error after removing the 2 outliers with an RMS error of0.88 Hz.

The estimation error can be further reduced by including the frequencytranslation and decimation step in (140) and (141) followed by theapplication of the AKF on the resulting signal. FIG. 51 plots themagnitude of the frequency estimation error for the case of SNR=−10 dBafter processing 280 samples corresponding to γ_(c) equal to 14.5 dB.The signal frequency selected randomly was −42.6 Hz for this case. Thereis no outlier in this case. The RMS frequency estimation error for thiscase is equal to 0.74 Hz. Examination of FIG. 51 shows that unlike thecase of Gaussian distribution, some of the errors are much higher than3×0.74=2.2 Hz. For a Gaussian distribution with RMS value of 0.74 Hz,the probability of having a value with its magnitude higher than 2.2 Hzis 0.0026 and thus among 100 simulation runs the expected number of runswith the magnitude of the estimation error exceeding 2.2 Hz is 0.26. Ifthe simulation runs with the magnitude of the estimation error exceeding7 times the C-R bound of 0.11 Hz is excluded from the statistics thenthe resulting RMS error is only 0.12 Hz for the 93 runs. FIG. 52 showsthe magnitude of the estimation error for those simulation runs forwhich the estimation error is less than 7 times the C-R bound off_(crb)=0.11. Thus, except for a few simulation runs, the RMS estimationerror is very close to the C-R bound. With the increase in the number ofsamples N_(s) the number of simulation runs for which the error exceeds7f_(crb) is expected to reduce.

Applications

A frequency estimation system is one of the most important subsystems ofcommunication, navigation, radar and various other engineering systems.In some cases, the efficient and precise estimation may be the criticalcomponent in the system design and may significantly limit theperformance of these systems with respect to various metrics ofperformance. For example, in application to high dynamic GPS receivers,the ability to acquire and track the GPS carrier signal under dynamicconditions limits the performance and applicability of these receiversto various important applications. The present disclosure candrastically increase the performance capability of such systems over theexisting ones.

In terms of communication systems, precise frequency and phase of thecarrier are important in communication systems involving coherentmodulation techniques such as MQAM and MPSK. In traditionalcommunication applications, either the techniques that use square lawloops or Costas type loops are used to derive the carrier frequency andphase from the modulated signal or a pilot signal is used which istracked. The use of square law loops or Costas type loops results insignificant loss in terms of phase noise of the reference carrier andphase ambiguity problems wherein there is phase ambiguity in the carrierphase equal to integer multiple of 2π/M for an MPSK signal. The use ofpilot carrier results in a loss of signal power because a significantpart of available power is used up in the pilot. The ability to providefast and accurate frequency and phase estimates at very low SNRs canreduce the loss due to pilot carrier to an insignificant value. Forexample with the pilot signal 10 dB below the modulated signal power,the loss is only 0.45 dB, operation at 14 dB below results in only 0.17dB loss compared to a loss of 1-3 dB in a traditional system. Thesenumbers are based on the assumption that the carrier frequencyuncertainty is nearly equal to the modulation signal bandwidth.

More recently precise and fast frequency acquisition and tracking havebecome very important with the evolution of the OFDM (OrthogonalFrequency Division Multiplexing) in mobile communication systems. TheOFDM modulation scheme offers several advantages. For example, it hasreduced problems of inter-symbol interference (ISI) caused by multipathpropagation. It has superior performance in selective fadingenvironments. Due to these advantages, OFDM has become part of variousimportant standards such as WiMax. However, because the OFDM is based onthe orthogonality among various subcarrier signals, it is very importantthat this orthogonality be maintained when these subcarriers arereceived at the receiver. However, the mobile wireless channelsintroduce frequency offsets which cause the disruption of theorthogonality among the subcarriers resulting in mutual interferenceamong the various subcarriers. Therefore, it is very important toprecisely estimate such frequency offsets and correct them to avoid theproblem of intercarrier interference. The offsets may be functions oftime and may be different for different subcarriers. Thus it isnecessary that precise estimation of the frequency offset be made withminimum requirements on the estimation time and SNR which is alsolimited in systems involving error correction coding techniques. Thedisclosure can virtually eliminate the need for pilot carriers that areprovided for the purpose of estimating the carrier offsets resulting inincreased efficiency and yet are capable of providing precise estimatesfor these offsets.

In various radar systems, for example in Doppler and FM chirp radars,the ability to provide accurate acquisition and tracking in very low SNRconditions can result in a drastic reduction in transmit power resultingin “quiet radar.” In certain embodiments, the systems and methodsdisclosed herein may be applied to various multi-tracking applicationsincluding, but not limited to RADAR multi-target tracking, SONARmulti-target tracking, LIDAR multi-target tracking, or combinationsthereof. The above only provides a few examples where the disclosure canbe exploited. The disclosure can be exploited in diverse applicationsincluding satellite communication, terrestrial wireless communication,digital TV, radars, broadcasting, radio astronomy, aeronautical andspace systems, structural vibrations, seismology, generalinstrumentation, etc.

In general, it will be apparent that the embodiments described hereinmay be implemented in many different embodiments of software, firmware,and/or hardware, for example, based on Field Programmable Gate Array(FPGA) chips or implemented in Application-Specific Integrated Circuits(ASICS). The software and firmware code may be executed by a computer orcomputing device comprising a processor (e.g., a DSP or any othersimilar processing circuit) including, for example, the computing devicedescribed below. The processor may be in communication with memory oranother computer-readable medium comprising the software code. Thesoftware code or specialized control hardware that may be used toimplement embodiments is not limiting. For example, embodimentsdescribed herein may be implemented in computer software using anysuitable computer software language type, using, for example,conventional or object-oriented techniques. Such software may be storedon any type of suitable computer-readable medium or media, such as, forexample, a magnetic or optical storage medium. According to variousembodiments, the software may be firmware stored at an EEPROM and/orother non-volatile memory associated with a DSP or other similarprocessing circuit. The operation and behavior of the embodiments may bedescribed without specific reference to specific software code orspecialized hardware components. The absence of such specific referencesis feasible, because it is clearly understood that artisans of ordinaryskill would be able to design software and control hardware to implementthe embodiments based on the present description with no more thanreasonable effort and without undue experimentation.

FIG. 53 shows an example of a computing device 5300 according to oneembodiment. For the sake of clarity, the computing device 5300 isillustrated and described here in the context of a single computingdevice. However, it is to be appreciated and understood that any numberof suitably configured computing devices can be used to implement adescribed embodiment. For example, in at least some implementations,multiple communicatively linked computing devices may be used. One ormore of these devices can be communicatively linked in any suitable waysuch as via one or more networks. One or more networks can include,without limitation: the Internet, one or more local area networks(LANs), one or more wide area networks (WANs) or any combinationthereof.

In the example of FIG. 53, the computing device 5300 comprises one ormore processor circuits or processing units 5302, one or more memoryand/or storage circuit component(s) 5304 and one or more input/output(I/O) devices 5306. Additionally, the computing device 5300 comprises abus 5308 that allows the various circuit components and devices tocommunicate with one another. The bus 5308 represents one or more of anyof several types of bus structures, including a memory bus or memorycontroller, a peripheral bus, an accelerated graphics port, and aprocessor or local bus using any of a variety of bus architectures. Thebus 5308 may comprise wired and/or wireless buses.

The processing unit 5302 may be responsible for executing varioussoftware programs such as system programs, application programs, and/orprogram modules/blocks to provide computing and processing operationsfor the computing device 5300. The processing unit 5302 may beresponsible for performing various voice and data communicationsoperations for the computing device 5300 such as transmitting andreceiving voice and data information over one or more wired or wirelesscommunications channels. Although the processing unit 5302 of thecomputing device 5300 is shown in the context of a single processorarchitecture, it may be appreciated that the computing device 5300 mayuse any suitable processor architecture and/or any suitable number ofprocessors in accordance with the described embodiments. In oneembodiment, the processing unit 5302 may be implemented using a singleintegrated processor. The processing unit 5302 may be implemented as ahost central processing unit (CPU) using any suitable processor circuitor logic device (circuit), such as a general purpose processor. Theprocessing unit 5302 also may be implemented as a chip multiprocessor(CMP), dedicated processor, embedded processor, media processor,input/output (I/O) processor, co-processor, microprocessor, controller,microcontroller, application-specific integrated circuit (ASIC), fieldprogrammable gate array (FPGA), programmable logic device (PLD), orother processing device in accordance with the described embodiments.

As shown, the processing unit 5302 may be coupled to the memory and/orstorage component(s) 5304 through the bus 5308. The bus 5308 maycomprise any suitable interface and/or bus architecture for allowing theprocessing unit 5302 to access the memory and/or storage component(s)5304. Although the memory and/or storage component(s) 5304 may be shownas being separate from the processing unit 5302 for purposes ofillustration, it is worthy to note that in various embodiments someportion or the entire memory and/or storage component(s) 5304 may beincluded on the same integrated circuit as the processing unit 5302.Alternatively, some portion or the entire memory and/or storagecomponent(s) 5304 may be disposed on an integrated circuit or othermedium (e.g., hard disk drive) external to the integrated circuit of theprocessing unit 5302. In various embodiments, the computing device 5300may comprise an expansion slot to support a multimedia and/or memorycard, for example.

The memory and/or storage component(s) 5304 represent one or morecomputer-readable media. The memory and/or storage component(s) 5304 maybe implemented using any computer-readable media capable of storing datasuch as volatile or non-volatile memory, removable or non-removablememory, erasable or non-erasable memory, writeable or re-writeablememory, and so forth. The memory and/or storage component(s) 5304 maycomprise volatile media (e.g., random access memory (RAM)) and/ornon-volatile media (e.g., read only memory (ROM), Flash memory, opticaldisks, magnetic disks and the like). The memory and/or storagecomponent(s) 5304 may comprise fixed media (e.g., RAM, ROM, a fixed harddrive, etc.) as well as removable media (e.g., a Flash memory drive, aremovable hard drive, an optical disk). Examples of computer-readablestorage media may include, without limitation, RAM, dynamic RAM (DRAM),Double-Data-Rate DRAM (DDRAM), synchronous DRAM (SDRAM), static RAM(SRAM), read-only memory (ROM), programmable ROM (PROM), erasableprogrammable ROM (EPROM), electrically erasable programmable ROM(EEPROM), flash memory (e.g., NOR or NAND flash memory), contentaddressable memory (CAM), polymer memory (e.g., ferroelectric polymermemory), phase-change memory, ovonic memory, ferroelectric memory,silicon-oxide-nitride-oxide-silicon (SONOS) memory, magnetic or opticalcards, or any other type of media suitable for storing information.

The one or more I/O devices 5306 allow a user to enter commands andinformation to the computing device 5300, and also allow information tobe presented to the user and/or other components or devices. Examples ofinput devices include data ports, analog to digital converters (ADCs),digital to analog converters (DACs), a keyboard, a cursor control device(e.g., a mouse), a microphone, a scanner and the like. Examples ofoutput devices include data ports, ADCs, DACs, a display device (e.g., amonitor or projector, speakers, a printer, a network card). Thecomputing device 5300 may comprise an alphanumeric keypad coupled to theprocessing unit 5302. The keypad may comprise, for example, a QWERTY keylayout and an integrated number dial pad. The computing device 5300 maycomprise a display coupled to the processing unit 5302. The display maycomprise any suitable visual interface for displaying content to a userof the computing device 5300. In one embodiment, for example, thedisplay may be implemented by a liquid crystal display (LCD) such as atouch-sensitive color (e.g., 76-bit color) thin-film transistor (TFT)LCD screen. The touch-sensitive LCD may be used with a stylus and/or ahandwriting recognizer program.

The processing unit 5302 may be arranged to provide processing orcomputing resources to the computing device 5300. For example, theprocessing unit 5302 may be responsible for executing various softwareprograms including system programs such as operating system (OS) andapplication programs. System programs generally may assist in therunning of the computing device 5300 and may be directly responsible forcontrolling, integrating, and managing the individual hardwarecomponents of the computer system. The OS may be implemented, forexample, as a Microsoft® Windows OS, Symbian OS™, Embedix OS, Linux OS,Binary Run-time Environment for Wireless (BREW) OS, JavaOS, or othersuitable OS in accordance with the described embodiments. The computingdevice 5300 may comprise other system programs such as device drivers,programming tools, utility programs, software libraries, applicationprogramming interfaces (APIs), and so forth.

In various embodiments disclosed herein, a single component may bereplaced by multiple components, and multiple components may be replacedby a single component to perform a given function or functions. Exceptwhere such substitution would not be operative, such substitution iswithin the intended scope of the embodiments.

While various embodiments have been described herein, it should beapparent that various modifications, alterations, and adaptations tothose embodiments may occur to persons skilled in the art withattainment of at least some of the advantages. The disclosed embodimentsare therefore intended to include all such modifications, alterations,and adaptations without departing from the scope of the embodiments asset forth herein.

Embodiments may be provided as a computer program product including anon-transitory machine-readable storage medium having stored thereoninstructions (in compressed or uncompressed form) that may be used toprogram a computer (or other electronic device) to perform processes ormethods described herein. The machine-readable storage medium mayinclude, but is not limited to, hard drives, floppy diskettes, opticaldisks, CD-ROMs, DVDs, read-only memories (ROMs), random access memories(RAMs), EPROMs, EEPROMs, flash memory, magnetic or optical cards,solid-state memory devices, or other types of media/machine-readablemedium suitable for storing electronic instructions. Further,embodiments may also be provided as a computer program product includinga transitory machine-readable signal (in compressed or uncompressedform). Examples of machine-readable signals, whether modulated using acarrier or not include, but are not limited to, signals that a computersystem or machine hosting or running a computer program can beconfigured to access, including signals downloaded through the Internetor other networks. For example, the distribution of software may be anInternet download.

Although embodiments have been described in language specific tostructural features and/or methodological acts, it is to be understoodthat the disclosure is not necessarily limited to the specific featuresor acts described. Rather, the specific features and acts are disclosedas illustrative forms of implementing the embodiments. Conditionallanguage, such as, among others, “can,” “could,” “might,” or “may,”unless specifically stated otherwise, or otherwise understood within thecontext as used, is generally intended to convey that certainembodiments could include, while other embodiments do not include,certain features, elements, and/or steps. Thus, such conditionallanguage is not generally intended to imply that features, elements,and/or steps are in any way required for one or more embodiments or thatone or more embodiments necessarily include logic for deciding, with orwithout user input or prompting, whether these features, elements,and/or steps are included or are to be performed in any particularembodiment.

That which is claimed:
 1. A method, comprising: deconstructing, by oneor more processors, a received input signal comprising a plurality ofconstituent signals, into a first plurality of component signals,wherein each of the first plurality of component signals includes acorresponding one or more second component signals; iterativelyperforming, by the one or more processors, a first adaptive linearestimation on each of the one or more second component signals togenerate a state vector corresponding to each of the one or more secondcomponent signals; iteratively performing, by the one or moreprocessors, one or more second adaptive linear estimations on each ofthe state vectors to generate a plurality of parameter vectorscorresponding to each of the state vectors; iteratively estimating, bythe one or more processors, based at least in part on the each of theplurality of parameter vector and its corresponding state vector, one ormore component parameters associated with the corresponding secondcomponent signal; combining, by the one or more processors, theestimated one or more component parameters to generate signal parametersassociated with the constituent signals; calculating, by the one or moreprocessors, respective correlations between the input signal andreference signals generated based at least in part on the signalparameters; and determining, by the one or more processors, that therespective correlations between the input signal and reference signalsmeet a threshold condition.
 2. The method of claim 1, wherein the inputsignal is a complex baseband signal that is a down conversion of andtime sampling of a radio frequency (RF) signal.
 3. The method of claim1, wherein the input signal is an Orthogonal Frequency Division MultipleAccess (OFDM) baseband signal comprising of a plurality of subcarriersignals and wherein estimating the signal parameters comprisesestimating at least one of: (i) a frequency offsets of the subcarriersignals; (ii) a phase of the subcarrier signals; or (iii) amplitudefades of the subcarrier signals introduced by the frequency selectivefading channel during the transmission of the OFDM baseband signal. 4.The method of claim 1, wherein performing the first adaptive linearestimation on each of the second component signals comprises using afirst adaptive Kalman filter and performing the one or more secondadaptive linear estimations on each of the state vectors comprises usinga second set of adaptive Kalman filters.
 5. The method of claim 1,wherein deconstructing the received input signal into the firstplurality of component signals comprises applying the input signal to atleast one of a bank of fast Fourier transform (FFT) filters or a bank ofpolyphase filters.
 6. The method of claim 1, wherein iterativelyestimating the one or more component parameters further comprises:generating a reference signal based at least in part on the one or morecomponent parameters corresponding to a second component of a firstcomponent signal of the plurality of component signals; calculating acorrelation between the reference signal and the first component signal;generating a new estimate of the one or more component parameters untilthe correlation is below a threshold.
 7. The method of claim 6, whereinthe threshold is based at least in part on the one or more componentparameters corresponding the second component signal.
 8. The method ofclaim 6, wherein the correlation is inversely related to an estimationerror between the first reference signal and a second component of thefirst component signal.
 9. The method of claim 6, wherein thecorrelation is performed between a reference signal based at least inpart on the one or more signal parameters corresponding to a firstconstituent signal and the input signal.
 10. The method of claim 1,further comprising generating a new estimates of the one or more signalparameters using a plurality of second stage adaptive Kalman filters.11. The method of claim 1, wherein the respective sets of firstcomponent signals include sinusoidal signals and the one or moreestimated signal parameters comprise at least one of: (i) one or moreamplitudes of the input signal; (ii) one or more phases of the inputsignal; (iii) one or more frequencies of the input signal; (iv) one ormore first derivatives of the frequencies of the input signal; (v) oneor more higher order derivatives of the frequencies of the input signal;or (vi) one or more slew rates of the input signal.
 12. The method ofclaim 1, wherein the input signal comprises a baseband version of areceived frequency chirp Radar signal.
 13. The method of claim 1,wherein combining the one or more estimated component parameters togenerate signal parameters associated with the input signal comprisesapplying the one or more estimated component parameters to a combiner orcollator.
 14. A system, comprising: a memory configured to store timeseries samples of an input signal comprising a plurality of constituentsignals; one or more filters for deconstructing the time series of inputsignal to a time series of a plurality of first component signals,wherein each of the plurality of first component signals comprise one ormore second component signals; a plurality of multi-frequency adaptiveKalman filters configured to receive the time series samples of thefirst component signals and generate at least one parameter estimatecorresponding to each of the one or more second component signals of thetime series associated with the first component signal; one or morecollators for combining the parameter estimates corresponding to each ofthe plurality of second component signals to generate estimates ofsignal parameters associated with the input signal; one or morecorrelation calculators for calculating the correlation between the timeseries samples of the input signal and a reference signal generatedbased at least in part on the signal parameter estimates; a plurality ofcorrelation threshold detectors, each correlation threshold detectorconfigured to compare a magnitude of each of the correlation from theone or more correlation calculators with corresponding respectivethreshold levels; an algorithm controller configured to generate atleast one control signal based at least in part on at least one of theplurality of correlation signals and provide the at least one controlsignal one of the (i) memory, or (ii) the plurality of multi-frequencyadaptive Kalman filter.
 15. The system of claim 14, wherein each of theplurality of multi-frequency adaptive Kalman filters comprise: a statevector Kalman filter configured to receive the time series samples ofthe constituent signals and generate a state vector time seriesassociated with the constituent signals; a vector splitter configured toreceive the state vector time series and generate a plurality of statevector time series subsets; a plurality of parameter vector Kalmanfilters, each parameter vector Kalman filter corresponding to arespective state vector time series subset and configured to generate aparameter vector time series based at least in part on the correspondingstate vector time series subset; a signal parameter estimator configuredto receive the plurality of parameter vector time series and the statevector time series to generate at least one estimated signal parametertime series associated with the input signal.
 16. The system of 15,wherein the signal parameter estimates are a first signal parameterestimates and further comprising a plurality of second stage adaptiveKalman filters, wherein each of the second stage adaptive Kalman filtersis configured to receive the first signal parameter estimates andgenerate a second signal parameter estimates with a reduced estimationerror.
 17. The system of claim 16, wherein each of the second stageadaptive Kalman filters comprise: a first Kalman filter configured toreceive a corresponding respective signal parameter estimate time seriesassociated with the input signal and generate a second state vector timeseries; a second Kalman filter configured to receive the second statevector time series and generate a second parameter vector time seriesbased at least in part on the received second state vector time series;a signal parameter estimator configured to receive the second parametervector time series and the second state vector time series and generateat least one estimated signal parameter associated with the inputsignal.
 18. The system of claim 14, wherein the input signal is anOrthogonal Frequency Division Multiple Access (OFDM) baseband signalcomprising of a plurality of subcarrier signals and wherein estimatingthe signal parameters comprises estimating at least one of: (i) afrequency offsets of the subcarrier signals; (ii) a phase of thesubcarrier signals; or (iii) amplitude fades of the subcarrier signalsintroduced by the frequency selective fading channel during thetransmission of the OFDM baseband signal.
 19. The system of claim 14,wherein each the threshold levels are based at least in part onestimates of the signal parameters associated with the plurality ofcomponent signals.
 20. The system of claim 14, wherein the correspondingrespective correlation signal is indicative of a convergence of theestimate of the signal parameters associated with the plurality ofconstituent signals.
 21. The system of claim 14, wherein the constituentsignals comprise sinusoidal signals and the one or more estimated signalparameters comprise at least one of: (i) one or more amplitudes of theinput signal; (ii) one or more phases of the input signal; (iii) one ormore frequencies of the input signal; (iv) one or more first derivativesof the frequencies of the input signal; (v) one or more higher orderderivatives of the frequencies of the input signal; or (vi) one or moreslew rates of the input signal.
 22. At least one non-transitorycomputer-readable medium comprising computer- executable instructionsthat, when executed by one or more processors, executes a methodcomprising: deconstructing a received input signal, comprising one ormore constituent signals, into one or more component signals; performinga first adaptive linear estimation on each of the component signals togenerate at least one state vector corresponding to each of thecomponent signals; separating elements of each of the at least one statevectors into one or more corresponding respective subsets; performingone or more second adaptive linear estimations on each of the subsets togenerate at least one parameter vector corresponding to each subset;estimating, based at least in part on the at least one parameter vectorand the at least one state vector, one or more signals parametersassociated with the constituent signal; calculating respectivecorrelations between the input signal and reference signals generatedbased at least in part on the one or more signal parameters; anddetermining that the respective correlations between the input signaland reference signals meet a threshold condition.
 23. The at least onenon-transitory computer-readable medium of claim 22, wherein the one ormore signal parameters comprise at least one of: (i) one or moreamplitudes of the input signal; (ii) one or more phases of the inputsignal; (iii) one or more frequencies of the input signal; (iv) one ormore first derivatives of the frequencies of the input signal; (v) oneor more higher order derivatives of the frequencies of the input signal;or (vi) one or more slew rates of the input signal.
 24. The at least onenon-transitory computer-readable medium of claim 22, further comprisingcombining the one or more signal parameters corresponding to thecomponent signals to generate the signal parameters associated with theconstituent signals.
 25. The at least one non-transitorycomputer-readable medium of claim 24, wherein combining the one or moresignal parameters corresponding to the component signals comprisesapplying the signal parameters corresponding to each of the componentsignals to a combiner or collator.